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Chapter 1 


INTRODUCTION 

The finite element method is one of the most significant develop- 
ments for solving problems of continuum mechanics. It was first 
applied by Turner et al. [l]* in 1956 for the analysis of complex 
aerospace structures. With increasing availability of digital 
computers, the method has become widespread and well recognized as 
applicable to a variety of continuum problems. Applications of the 
method to thermal problems were introduced in the middle of 1960 T s 
for the solution of steady-state conduction heat transfer [2]. 
Thereafter, extensions of the method were made to both transient and 
nonlinear analyses where nonlinearities may arise from temperature 
dependent material properties and nonlinear boundary conditions. 
Important publications of finite element heat transfer analysis 
appear in references [3-12], With these developments and consider- 
able effort contributed during the past decade, the method has 
gradually increased in thermal analysis capability and become a 
practical technique for analyzing realistic thermal problems. 


*The numbers in brackets indicate references. 



1.1 Current Status of Thermal-Structural Analysis 


Thermal stresses induced by aerodynamic heating on advanced 
space transportation vehicles are an important concern in structural 
design. Nonuniform heating may have a significant effect on the 
performance of the structures and efficient techniques for determining 
thermal stresses are required. Frequently, the thermal analysis of 
the structure is performed by the finite difference method. 
Production-type finite difference programs such as MITAS and SINDA 
have demonstrated excellent capabilities for analyzing complex 
structures [13] . In structural analysis, however, the finite element 
method is favorable due to better capabilities in modeling complex 
structural geometries and handling various types of boundary condi- 
tions. To perform coupled thermal-structural analysis with efficiency, 
a computer program which includes both thermal and structural analysis 
codes is preferred, and a single numerical method is desirable to 
eliminate the tedious and perhaps expensive task of transferring 
data between different analytical models. 

Currently, the capabilities and efficiency of the finite element 
method i^ analyzing typical heat transfer problems such as combined 
conduction-forced convection is about the same as using the finite 
difference method [14]. With the wide acceptance of the finite 
element method in structures and its rapid growth in thermal analysis, 
it is particularly well-suited for coupled thermal-structural 
analysis. At present, several finite element programs which include 
both thermal and structural analysis capabilities exist; e.g. 

NASTRAN, ANSYS, ADINA and SPAR are widely used. These programs use 
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a common data base for transferring temperatures computed from a 
thermal analysis processor to a structural analysis processor for 
determining displacements and stresses. With the use of a common 
finite element discretization, a significant reduction of effort in 
preparing data is achieved and errors that may occur by manually 
transferring data between analyses is eliminated. 

1.2 Needs for Improving Finite Element Methodology 

Although the finite element method offers high potential for 
coupled thermal-structural analysis, further improvements of the 
method are needed. Quite often, the finite element thermal model 
requires a finer discretization than the structural model to compute 
the temperature distribution accurately. Detailed temperature 
distributions are necessary for the structural analysis to predict 
thermal stress distributions including critical stress locations 
accurately. Improvement of thermal finite elements is, therefore, 
required so that a common discretization between the two analytical 
models can be maintained. 

Another need for improving the method includes a capability of 
the thermal analysis to produce thermal loads required for the 
structural analysis directly. At present, typical thermal-structural 
finite element programs only transfer nodal temperatures computed 
from the thermal analysis to the structural analysis. These nodal 
temperatures are generally inadequate because additional information, 
such as element temperature distributions and temperature gradients, 
may be required to compute thermal stress distributions correctly. 



4 


These needs are important in improvement of finite element 
coupled thermal-structural analysis capability. The use of improved 
thermal finite elements can reduce model size and computational 
costs especially for analysis of complex aerospace vehicle structures. 
Improved thermal elements will also have a direct effect in increasing 
the structural analysis accuracy through improving the accuracy of 
thermal loads. 

To meet these requirements for improved thermal-structural 
analysis and to demonstrate benefits that can be achieved, this 
dissertation will develop an approach called integrated finite 
element thermal-structural analysis. First, basic concepts of the 
integrated finite element thermal-structural formulation are intro- 
duced in Chapter 2. Finite elements which provide exact solutions 
to one-dimensional linear steady-state thermal-structural problems 
are developed in Chapter 3. Chapter 4 demonstrates the use of these 
finite elements for linear transient analysis. Next, in Chapter 5 
a generalized approach for improved finite elements is established 
and its efficiency is demonstrated through thermal-structural 
analysis with radiation heat transfer. Finally, in Chapter 6 
extension of the approach to two dimensions is made with a new two- 
dimensional finite element. In each chapter, benefits of utilizing 
the improved finite elements are demonstrated by both academic and 
realistic thermal-structural problems. 

Throughout the development of the improved finite elements, 
detailed analytical and finite element formulations are presented. 

Such details are provided in the form of equations, finite element 
matrices in tables and computer subroutines in appendices. 



Chapter 2 


AN INTEGRATED THERMAL- STRUCTURAL FINITE 
ELEMENT FORMULATION 

2,1 Basic Concepts 

Before applying the finite element method to thermal-structural 
analysis, it is appropriate to establish basic concepts and procedures 
of the method. Briefly described, the finite element method is a 
numerical analysis technique for obtaining approximate solutions to 
problems by idealizing the continuum model as a finite number of 

discrete regions called elements. These elements are connected at 

/ 

points called nodes where normally the dependent variables such as 
temperature and displacements are determined. Numerical computations 
for each individual element generate element matrices which are then 
assembled to form a set of linear algebraic equations (for 
steady state problems) to represent the entire problem. These 
algebraic equations are solved simultaneously for the unknown 
dependent variables. Usually the more elements used, the greater 
the accuracy of the results. Accuracy, however, can be affected by 
factors such as the type of element selected to represent the con- 
tinuum, and the sophistication of element interpolation functions. 
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2.2 Element Interpolation Functions 

The first step after replacing the continuum model by a 
discrete number of finite elements is to determine a functional 
relationship between the dependent variable within the element and 
the nodal variables* The function that represents the variation of 
a dependent variable is called the interpolation function. In thermal 
analysis, the element temperature T(x,y,z,t) are generally expressed 
in the form 

T(x,y,z,t) = [N T (x,y,z)J (T(t)} (2.1) 

where [N T (x,y,z)J denotes a row matrix of the element temperature 
interpolation functions, and { T ( t ) } denotes a vector of nodal 
temperatures. Similarly, in a structural analysis, the element 
displacements, {6} , are expressed as, 

{ 6 (x,y,z, t) } = [N g (x,y,z)] (5 (t) } (2.2) 

where [N^(x,y,z)] denotes a matrix of structural displacement inter- 
polation functions, and (6(t)} denotes a vector of nodal 
displacements. 

Usually, polynomials are selected as element interpolation 
functions and the degree of the polynomial chosen depends on the 
number of nodes assigned to the element. Regardless of the algebraic 
form, these interpolation functions have a value of unity at the node 
to which it pertains and a value of zero at other nodes. For example, 
linear temperature variation for a two-node one-dimensional rod 
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element with nodal temperatures and T ^ at x =0 (node 1) and 

x =L (node 2), respectively, can be written in the form 


T(x,t) - Ll - - 



T-^t) 


< > 


T 2 (t) 

/ 


By comparing this equation with the general form of the element 
temperature variation, Eq. (2.1), the element interpolation functions 
are 


N^x) =■ 1 - and N 2 (x) - ^ 


These element interpolation functions, therefore, have the properties 
of N. * 1 at node i and N. * 0 at the other node. 

i i 


2.3 Finite Element Thermal Analysis 

Once the type of elements and their interpolation functions have 
been selected, the matrix equations expressing the properties of the 
individual element are evaluated. In thermal analysis, the method of 
weighted residuals [15] is frequently employed starting from the 
governing differential equations. For condution heat transfer in a 
three-dimensional anisotropic solid ft bounded by surface T 
(Fig. 1), an energy balance on a small element is given by, 





+ Q(x,y,z,t) 


9T(x,y,z,t) 
3 1 


(2.3) 


where q , q , q are components of the heat flow race per unit area, 
x y z 

Q is the internal heat generation rate per unit volume, p is the 
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density, and c is the specific heat. Using Fourier’s Law, the 
components of heat flow rate for an anisotropic medium can be written 
in the matrix form. 




— 


— 


r 

q x 


k n 

k 12 

k 13 


9T/3x 

q y 

> m - 

k 

21 

k 22 

k 23 

< 

3T/3y > 

q z 


k 31 

k 32 

k 33 


3T/3z 

^ J 





_ 


^ J 


where is the symmetric conductivity tensor. Figure 1 shows 

several types of boundary conditions frequently encountered in the 
analysis. These boundary conditions are (1) specified surface 
temperatures, (2) surface heating, (3) surface convection, and 
(4) surface radiation: 


T - T S 


on 


(2.5a) 


qn + qn + qn = - q 
xx y y z z ^s 


qn + qn 4* q n = h(T c 
xx y y z z s 


qn + qn + qn = oz T 
xx y y z z s 


on S 2 
TJ on S 3 

aq r on S 4 


(2.5b) 


(2.5c) 


(2.5d) 


where T s is the specified surface temperature; n x , n^, n^ are the 
direction cosines of the outward normal to the surface, q s is the 
surface heating rate unit area, h is the convection coefficient, 

is the convective medium temperature, a is the Stef an-Boltzmann 
constant, e is the surface emissivity, a is the surface absorp- 
tivity, and .q r is the incident radiant heat flow rate per unit 


area. 
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To apply the finite element technique, the domain SI is first 
discretized into a number of elements. For an element with r 
nodes, the element temperature, Eq. (2.1), can be written in the 
form 


r 

T(x,y,z,t) = Z N. (x,y,z) T (t) (2.7a) 

i=l 

and the temperature gradients within each element are 


3T(x,y,z,t) 

3x 


r 

Z 

i=l 


3N i (x,y,z) T ± (t) 
3x 


(2.7b) 


3T(x,y,z,t) 

3y 


r 

Z 

i=l 


3N i (x,y,z) 

3Y 


T ± (t) 


(2.7c) 


3T(x,y,z,t) 

dz 


J 3N i (x,y,z) T ± (t) 
i=l 3z 


(2. 7d) 


These element temperature gradients can be written in the matrix 
form, 


_3T (x,y,z,t) 
3x 


3T (x,y,z,t) I 

3 y f 


3T (x,y,z,t) 

8 Z 

^ ) 


[B(x,y,z)] (T(t) } 


where 


[B(x,y,z)] 


is the temperature-gradient interpolation 


( 2 . 8 ) 


matrix: 
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[B(x,y,z)] 


9N X 3N2 ^ 3N r 

3x 3x 3x 


SNj^ 9N 2 9N r 

9y 9y 9y 


9N 1 3N 2 3N r 

3z 3z 3z 


(2.9) 


and, therefore, the components of heat flow rate, Eq. (2.4), become 



-[k] [B] {T } 


( 2 . 10 ) 


where [k] denotes the thermal conductivity matrix. 

In the derivation of the element equations, the method of 
weighted residuals is applied to the energy equation, Eq. (2.3), for 
each individual element (e) . This method requires 


J 


(e) 


9 q 3q 


9x 


3y 



J 9T s 

M at 2 


N.dft = 0 


(2.11) 


1 “ 1)2 ■•••IT 


After the integrations are performed on the first three terms by 

using Gauss’s Theorem, a surface integral of the heat flow across 

( s ) 

the element boundary, T , is introduced, and the above equations 
become 


( [ (q -n) N dT 

V J r (e) 







'i 

q x 


r 

3N ± 

3N. 

3N. 


n 

> dn J 


9x 

9y 

3z _ 

< 

q y 

q z 

^ J 



12 


r 


.(e) QN i d " + J, 


, . pc 4? N. d.Q = 0 
(e) at 1 


( 2 . 12 ) 


where q is the vector of conduction heat flux across the element 
boundary and n is a unit vector normal to the boundary. The 
boundary conditions as shown in Eqs. (2.5a -2.5d) are then imposed, 



(q • n) dr 


f q N. dr + 
I s i 




h(T -TJ N ± 


dr 


+ I (aeT 

s. 


aq f ) dr 




(e) L 


3N . 3N. 3N. 

i x x 

3x 3y 3z 


> dft 


,(e) 


QN ± dfi + 


(e) 


pc 4? N - 
3t x 


= 0 


(2.13) 


By substituting the vector of heat flow rate, Eq. (2.10), the above 
element equations finally result in the matrix form, 


[C] { T) + [[K c ] + [1^] + [Kj]{T} 

= (R c ) + { R q ) + { R q } + {R h } + {R r } 


(2.14) 


where [c] is the element capacitance matrix; [K^] , [K^] and 

[K^] are element conductance matrices corresponding to conduction, 
convection and radiation, respectively. These matrices are expressed 


as follows : 
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[c] 

= J n (e) pc {N T } W dn 

(2.15a) 


[K C J 

= | (e) [B t ] T [k][B T ] dfl 

(2.15b) 


'V 

S 3 

h {n t > |n t J dr 

(2.15c) 


[kJ {t} 

■J. 

b 4 

aeT 4 {N t > dr 

(2 . 15d) 

The right-hand side of 

the discretized equation 

(2.14) contains 

heat load vectors due to specified nodal temperatures 

, internal heat 

generation, 

specified surface heating, surface convection and surface 

radiation. 

These vectors are defined by 



u c ) 

= - 

(q-n) (N } dr 

S 1 

(2.16a) 


<y 

0 

[(e) Q { V d!! 

(2.16b) 


{ V 

= 

L % { V dr 
s 2 

(2.16c) 


{ V 

J 

r 

h t m { N } dr 
S 3 

(2 . 16d) 



- J 

aq { N } dr 

Q 1 

(2 . 16e) 


where q is the vector of conduction heat flux across boundary that 
is required to maintain the specified nodal temperatures. 
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2.4 Finite Element Structural Analysis 

In a finite element structural analysis, element matrices may 
be derived by the method of weighted residuals, or by a variational 
method such as the principle of minimum potential energy [ 17—19 3 . 

For simplicity in establishing these element matrices and understand- 
ing general derivations, the last approach is presented herein. 

The basic idea of this approach is to derive the static equilibrium 
equations and then include dynamic effects through the use of 
D'Alembert's principle. 

Consider an elastic body in a three-dimensional state of stress. 
The internal strain energy of an element (e) can be written in a 
form, 

u = \ £ (e) |e " £ 0 J {a} (2 * 17) 

C s ) 

where Q'’ is the element volume, {a} denotes a vector of stress 

components; [eJ and [e^J denote row matrices of total strain and 

initial strain components, respectively. Using the stress-strain 
relations, 

{a} - [D] { e - e Q } (2.18) 

where [D] is the elasticity matrix, the internal strain energy 
becomes 

D - 7 J ( S ) - e cJ [D] te - e 0 > d “ 
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or 


U-i 


(e) l £ J [D] {e} dft - J (e) |_ej [D] {e Q } 




+ I | Q (e) ^ [D] {e 0 } d “ 


( 2 . 19 ) 


For each element, the potential energy of the external forces 
may result from body forces and boundary surface tractions. The 
potential energy due to body forces can be written in a form. 


V B - - L.) W «> da 

M i 


( 2 . 20 ) 


where {f} denotes a vector of body force components. Similarly, 
the potential energy due to surface tractions is. 


v s - - L.) l*j dr 


( 2 . 21 ) 


where {g} denotes a vector of surface traction components, and 

( e ) 

T denotes the element boundary. The total element potential 

energy, tt^ , the sum of the internal strain energy and the potential 
energy of the external forces is, 

"e " J ( ( e ) l E J ( £ > dn ~ f (e ) l £ J (e 0 > dfi 
+ I f (e) l £ 0 J tD) (E 0 } d!i ' I (e) L«J {£) d!! 

db JwJ 

" [ (e) L«J U} dr ( 2 - 22 > 
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For a three-dimensional finite element with r nodes, the 
displacement field can be expressed as 


r -n 


C 

r 

u(x,y,z,t) 


I N (x,y,z) u. (t) 
i=l 

r 

v(x,y,z,t) 

► = ^ 

E N (x,y,z) v (t) > 
i=l L 

r 

w(x,y,z,t) 


E N. (x,y,z) w (t) 
i=l 


(2 . 23) 


where u, v, w are components of displacement in the three coordinate 
directions. The vector of strain components can be computed from 


{e} 


r "\ 



e 


au 

X 


ax 

£ 

y 


3v 

3y 

£ 


aw 

Z 


3z 

< 

> = < 


Y 


au | av 

xy 


ay ax 

Y 


av + aw 

yz 


3z 3y 

Y 


3u aw 

xz 

v. > 


3z 3x 

J 


(2 . 24) 


where [b ] is the strain-displacement interpolation matrix. By 
substituting the element displacement vector, Eq. (2.23), and the 
vector of strain components, Eq. (2.24) into Eq. (2.22), the total 
element potential energy is expressed in terms of the nodal displace- 
ment vector {6} as 
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* e - \ [<sj | (e) [B s ] T [D] [B s 3 dfi { 6 } - {6} f (e) [B s ] T CD] {e 0 } dft 

+ ~2 J" (e) t- e O-l ^ ^ J (e) dft 

- { 5 } | (e) [N s ] T {g} dr (2.25) 

The principle of minimum potential energy requires. 


3* e 

— J- = 0 

a{6} 

which yields the element equilibrium equations, 

[ K s ] {6} = {F c } + {F b } + {F s } + {F t } (2.26) 

where [ K g ] is the element stiffness matrix defined by 

[K s 3 » J (e) [B s ] T [D] [B s ] d^ (2.27a) 

£2 


The right hand side of the equilibrium equations contains force 
vectors due to concentrated forces, body forces, surface tractions 
and initial strain, respectively. The nodal force vectors due to 
body forces and surface tractions are 


<v - 


I a (.) [n sI t W 

J r (e) t N s] T <a> dr 


(2.27b) 


{ F_ } = 


(2.27c) 
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For initial strains from thermal effects, the corresponding nodal 
vector {F^} is due to the change of temperature from a reference 
temperature of the zero-stress state and may be written as 

{F T } = j (e) [B s ]T [d] {a} (T " T ref ) (2.27d) 

where {a} is a vector of thermal expansion coefficients, T is 
the element temperature distribution, and T re f is the reference 
temperature for zero stress. 

For elastic bodies subjected to dynamic loads, the effects of 
inertia and damping forces must be taken into account. Using 
D r Alembert's principle, the inertia force can be treated as a body 
force given by 

{f} = - p { 6 } (2.28a) 

where p is the mass per unit volume. By using element displacement 
variations, Eq. (2.23), this inertia force is expressed in terms of 
nodal displacements as 

(f) - - P [N g ] {6} (2.28b) 

Similarly, the damping force which is usually assumed to be propor- 
tional to the velocity can be expressed in the form, 

{ f ) = - V [ N s ] { 6 } (2.28c) 

where u is a damping coefficient. By substituting these inertia 
and damping forces, Eqs. (2.28b -2.28c), into Eq. (2.27b), the equi- 
valent nodal body forces shown in Eq. (2.27b) become 
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{F B> = 1(e) [N s ]T {f} " J Q(e) * [N s ] dfl {6} 

- J (e) [N s ] T p [N g ] dfi {£} (2.29) 

Finally, by using the static equilibrium equations, Eq. (2.26), with 
the above equivalent nodal body force, the basic equations of struc- 
tural dynamics can be written in the form, 

[M] {6} + [C s ] {6} + [K g ] {6} = {F c } + {F b } + {F g } + (F t > (2.30) 

where [m] and [C g ] are the element mass and damping matrices, 
respectively, and defined by 

[M] = J [N s ] T P [N s ] dfi (2.31a) 

Q 

[C s ] " [ (e) [N s ] T U [N s ] dfl (2.31b) 

In a general formulation of transient thermal-stress problem, 
the heat conduction equation (2.3) contains a mechanical coupling 
term in addition [16]. This coupling term represents the mechanical 
energy associated with deformation of the continuum and in some 
highly specialized problems (see Ref. 16) can affect the temperature 
solution. In most of engineering applications, fortunately, this term 
is insignificant and is usually disregarded in the heat conduction 
equation. This simplification permits transient thermal solutions 
and dynamic structural responses to be computed independently. 

For a structural analysis where the inertia and damping effects 
are negligible, the static structural response, Eq. (2.26), can be 



20 


computed at selected times corresponding to the transient thermal 
solutions. Such a sequence of computations, widely used in thermal- 
structural applications, is called a quasi-static analysis. Results 
of temperatures directly enter the structural analysis through the 
computation of the thermal nodal force vector, Eq. (2.27d). Tempera- 
tures also have an indirect effect on the analysis through the 
structural material properties, since the elasticity matrix [D] 
and the thermal expansion coefficient vector {a} are, in general, 
temperature dependent. Temperature dependent properties may result 
in a variation of the structural element stiffness matrix, Eq. (2.27a), 
throughout the transient response. 

2.5 Integrated Approach 

The representation of the element temperature distribution in 
the computation of structural nodal forces is an important step in 
the coupled thermal-structural finite element analysis. In typical 
production-type finite element programs, element nodal temperatures 
are the only information transferred from the thermal analysis to 
the structural analysis. This general procedure is shown schemati- 
cally in Fig. 2(a) and herein is called the conventional finite 
element approach. Since the conventional thermal analysis only 
provides nodal temperatures, an approximate temperature distribution 
is assumed in the structural analysis which results in a reduction 
in accuracy of displacements and thermal stresses. 

To improve the capabilities and efficiency of the finite 
element method, an approach called integrated thermal-structural 
analysis is developed as illustrated by Fig. 2(b). The goals of 
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INTEGRATED 



• THERMAL AND STRUCTURAL • IMPROVED THERMAL ELEMENTS 
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• NODAL TEMPERATURES 1 T I • COMPATIBLE THERMAL DATA TRANSFER 

TRANSFERRED AS REQUIRED BY STRUCTURAL ELEMENT 

• THERMAL FORCES BASED • THERMAL FORCES BASED ON ACTUAL 

ONLY ON NODAL IT} TEMPERATURE DISTRIBUTIONS 


(a) Conventional analysis (b) Integrated analysis 


Fig. 2. Conventional versus integrated thermal and structural analysis. 


22 


the integrated approach are to: (1) provide thermal elements which 

predict detailed temperature variations accurately, (2) maintain 
the same discretization for both thermal and structural models with 
fully compatible thermal and structural elements, and (3) provide 
accurate thermal loads to the structural analysis to improve the 
accuracy of displacements and stresses. 

These goals of the integrated approach require developing new 
thermal finite elements that can provide higher accuracy and effi- 
ciency than conventional finite elements. The basic restriction on 
these new thermal elements is the required compatability with the 
structural elements to preserve a common discretization. Detailed 
temperature distributions resulting from the improved thermal finite 
elements can provide accurate thermal loads required for the 
structural analysis by rigorously evaluating the thermal load 
integral, Eq. (2.27d). 



Chapter 3 


EXACT FINITE ELEMENTS FOR ONE-DIMENSIONAL 
LINEAR THERMAL-STRUCTURAL PROBLEMS 

In general, polynomials are selected as element interpolation 
functions to describe variations of the dependent variable within 
elements. In one-dimensional analysis, the simplest polynomial which 
provides a linear variation within an element is of the first order, 

4) = C X + C 2 x (3.1) 

where <f> denotes the dependent variable such as temperature or 
displacement; C^ and C 2 denote constants, and x is the coor- 
dinate of a point within the element. A finite element with two 
nodes is formulated by imposing the conditions at nodes, 

<Kx =0) = <f >* I <j> (x = L) = <t> 2 (3.2) 

where L is the element length; 4>-|_ and <f> 2 are nodal values at 
node 1 and 2, respectively. The dependent variable, therefore, can 
be written in terms of nodal values as 

* - d - f> +1 + <T> *2 

or in the matrix form. 
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I 
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= L N J {*} (3.3) 

where [nJ is the row matrix of element interpolation functions. 

The type of finite element where the dependent variable is assumed 
to vary linearly between the two element nodes is often used in one- 
dimensional problems and is called a conventional finite element 
herein, With the linear approximation, a large number of elements 
are required to represent a sharply varying dependent variable. In 
some special cases, however, conventional finite elements can provide 
exact solutions when the solutions to problems are in the form of a 
linear variation. For example, a linear temperature variation is the 
exact solution of one-dimensional steady-state heat conduction in a 
slab; therefore, the use of the conventional finite element leads to 
an exact solution. Further observation [20] has shown that, under 
some conditions, exact nodal values are obtained through the use of 
this element type. Temperatures for steady-state heat conduction 
with internal heat generation in a slab and deformations of a bar 
loaded by its own weight are examples of this case. In the past, 
the capability of conventional finite elements to provide exact 
solutions has been regarded as a property of the .particular equation 
being solved and not applicable to general problems. 
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In this chapter, finite elements that provide exact solutions 
to one-dimensional linear steady-state thermal-structural problems 
are given. The fundamental approach in developing exact finite 
elements is based on the use of exact solutions to one-dimensional 
problems governed by linear ordinary differential equations. A 
general formulation of the exact finite element is first derived and 
applications are made to various thermal-structural problems. 
Benefits of utilizing the exact finite elements are demonstrated 
by comparison with results from conventional finite elements and 
exact solutions. 


3.1 Exact Element Formulation 


In this section, a general derivation of exact finite elements 
is given. Exact finite elements for various thermal and structural 
problems are derived and described in detail in the subsequent 
sections. Consider an ordinary, linear, nonhomogeneous differential 
equation. 


n, .n-1. 

a d_i + a d — i + 

n , n n-1 , n-1 

dx dx 


a i H + V = r(x) 


(3.4) 


where x is the independent variable, 4>(x) is the dependent 
variable, a^, i=0, n are constant coefficients, and r(x) is 

the forcing function. A general solution to the above differential 
equation has the form 

n 

<Kx) = Z C. f . (x) + g (x) 
i=l 1 1 


(3.5) 
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where Ck are arbitrary constants, f^(x) are typical functions in 
the homogeneous solution and g(x) is a particular solution. For 
example, a typical one-dimensional steady-state thermal analysis is 
governed by second order differential equation of the form of 
Eq. (3.4) and has a general solution 

4>(x) = Cjl f*|_(x) + C 2 f 2 (x) + g(x) (3.6) 

By comparing this general solution with the solution in the form of 
polynomials used to describe a linear variation of dependent variable 
in the conventional finite element, Eq. (3.1), basic differences 
between these two solutions are noted: (1) the function f^(x) in 

the general solution to a given differential equation can be forms 
other than the polynomials, and (2) the general solution contains a 
particular solution g(x) which is known in general and depends on 
forcing function r(x) on the right hand side of the differential 
equation (3.4). 

3.1.1 Exact Element Interpolation Functions and Nodeless 
Parameters 

Once a general solution to a given differential is obtained, 
exact element interpolation functions can be derived. For a typical 
finite element with n degrees of freedom, n boundary conditions 
are required. With the general solution shown in equation (3.5), 
the required boundary conditions are 

(x i ) = <t> i i = 1, 2 n (3.7) 

where x. is the nodal coordinate and d> . is the element nodal 

i l 
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unknown at node i. After applying the boundary conditions, the 
exact element variation of 4>(x) has the form, 

n 

<J> (x) = G(x) + E N ± 4> i 
i-1 

where N^(x) is the element interpolation function corresponding 
to node i. The function G(x) is a known function associated with 
the particular solution. In general ,' this function can be expressed 
as a product of a spatial function Nq(x) and a scalar term <\>q 
which contains a physical forcing parameter such as body force, 
surface heating, etc.; 


G(x) = Nq (x) <J> 0 


and, therefore, the exact element <j>(x) variation becomes, 


n 

<f>(x) = Nq (x) <j> 0 + Z N i <t> i (3.8a) 

i=l 


or in the matrix form 


<Kx) - L N o N x N 2 .... N n J = 

Note that the element interpolation function N^(x^) has a value 
of unity at node i to satisfy the boundary conditions, Eq. (3.7), 
thus the spatial function N^(x) must vanish at nodes. Since the 


K 1 

>2 


\ - L n o n J 


^ o 


(3.8b) 
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term <(>q is a known quantity and neither relates to the element 
nodal coordinates nor is identified with the element nodes, it is 
called a nodeless parameter. Likewise, the corresponding spatial 
function Nq(x) is called a nodeless interpolation function. 
Comparison between element variations of a typical nodeless para- 
meter finite element, Eq. (3.8), and the conventional linear finite 
element, Eq. (3.3), is shown in Fig. 3. 


3.1.2 Exact Element Matrices 

After exact element interpolation functions are obtained, the 
corresponding element matrices can be formulated. For the governing 
ordinary differential equation, Eq. (3.4), typical element matrices 
can be derived (see section 2.2) and element equations can be written 
in the form. 



(3.9) 


where , i, j =0, n are typical terms in the element stiffness 

matrix; F^, i =0, n are typical terms in the element load vector, 
<t>^, i =1, n are the element nodal unknowns, and <p^ is the element 
nodeless parameter. Since the element nodeless parameter is known. 


the above element equations reduce to 
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X 


Fig. 3. Comparison of conventional and nodeless parameter 
elements. 
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(3.10) 


3.2 Exact Finite Elements in Thermal Problems 

In one-dimensional linear steady-state thermal problems, typical 
governing differential equations can be derived from a heat balance 
on a small segment in the form, 

d! (a 2 (x) + a l (x) £ + a 0 (x) T = r(x) (3 ' 1X) 

where T denotes the temperature, x denotes a typical one- 
dimensional space coordinate in cartesian, cylindrical or spherical 
coordinates; a^, i =0, 1,2 are variable coefficients, and r(x) 
is a function associated with a heat load for a given problem, A 
general solution to the above differential equation has the form, 

T (x) = C 1 f 1 (x) + C 2 f 2 (x) + g(x) (3.12) 

where f^(x) and f 2 (x) are -*-^ near ^y independent solutions of the 
homogeneous equation, and C 2 are constants of integration, 

and g(x) is a particular solution. Since the particular solution 
g(x) is known, the above general solution has two unknowns to be 
determined. A finite element with two nodes, therefore, can be 
formulated using the conditions. 
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T( X;l ) = T x (3.13a) 

T(x 2 ) = T 2 (3.13b) 

where x., i=l, 2 are nodal coordinates and T., i= 1,2 are the 
x* x* * 

nodal temperatures. Imposing these conditions on the general solu- 
tion yields two equations for evaluating and C^, 

T(x 1 ) = f^x^ + C 2 f 2 ( X;L ) + g( X;L ) 

X(x 2 ) = t 2 = c i f i^ x 2^ + C 2 f 2 < - x 2^ + 8 ^ x 2^ 

or in matrix form 


w 

f 2 ( x 1 ) 

< 

f •> 

C 1 

► = 1 

r 

T x - g(x x ) 

> 

fj_(x 2 ) 

f 2 ( x 2 ) 


C 2 


T 2 " g(x 2 ) 

^ J 


After and are determined and substituted into the general 

solution, Eq. (3.12), the exact element temperature variation can 
be written as 


T (x) = N q (x) T q + N x (x) T x + N 2 (x) T 2 


(3.14a) 


or in the matrix form. 


r 

Tn 


T(x) = LNq \ N 2 J 1 T 1 

T 


(3.14b) 


where Nq(x) is the nodeless interpolation function and Tq 


is the 
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nodeless parameter; N-^(x) and N^(x) are element interpolation 
functions corresponding to node 1 and 2, respectively. These element 
interpolation functions including the nodeless parameter are known 
functions defined by 


f,(x) f 9 (x 2 ) - f 7 0O f 9 (x) 

N (x) = — - - (3.15a) 

W 

f-^x^) f 2 (x) - f^x) f 2 (x 1 ) 

N 2 (x) = (3.15b) 

W 


N 0 (x) 


f,(x.) g(x_) - f,(x,) g(x.) 

T Q = g(x) + — — f. (x) 

W 1 


f 1 (x.) g (x- ) - f. (x ) g(x„) 

+ — f 2 (x) (3.15c) 

W 


where W = f^x^ f 2 (x 2 ) ” f i^ x 2^ f 2^ x l^ ' 


Using the exact element interpolation functions shown in 
Eq. (3.14), and the governing differential equation, Eq. (3.11), 
element matrices can be derived through the use of the method of 
weighted residuals; 


r d , dT. 
dx a 2 dx 


dT . 

a i 1 l " 

1 dx 


£ + a Q X " r ] 


N. dx = 0 
1 


i=0,l,2 


(3.16) 


Performing an integration by parts on the first term and substituting 
for element temperature in terms of the interpolation functions, 

Eq. (3.14), yields element equations in the form. 
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[[K c ] 


+ [ig + [k h ]] ^ t 


> = (Q r > + {Q} 


(3.17) 


where [K c ], [Ky] , and [KjJ are the element conductance matrices 
associated with the second, the first, and the zero-order derivative 
term on the left hand side of the governing differential equation 
(3.11), respectively; {Q c > is the element vector of conduction 
heat flux across element boundary, and {Q} is the element load 
vector from the heat load r(x") in the governing differential 
equation. These matrices are defined as follows: 



(3.18a) 


(3.18b) 


x„ 


'V - 


1 Q {N} [nJ dx 


(3.18c) 


{Q c > - < 


r dT 

I' a 2 dJ N 1 


(3 . 18d) 


I 
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(Q) 



(3 . 18e) 


Depending on the complexity of element interpolation functions, the 
element matrices may be evaluated in closed form or they may require 
numerical integration. However, after the element matrices are 
computed, typical element equations can be written in the form. 



(3.19) 


Since the nodeless parameter is known, the first equation is uncoupled 
from the nodal unknowns in the second and third equations. Thus, 
the exact element matrices have the same size as of the conventional 
linear finite element and element equations can be written as 



Note that, in general, the above conductance matrix is an 
asymmetric matrix. This asymmetry is caused by the conductance 
matrix [K ] shown in Eq. (3.18b) associated with the first-order 
derivative in the governing differential equation, Eq. (3.11). To 
obtain a symmetrical conductance matrix, the first-order derivative 
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is eliminated by casting the governing differential equation in self- 
adjoint form. 


di C p(x > + Q(x) T = R(x) 


(3.21) 


where 


P(x) = expcjfcaj^ + -^) /a 2 ] dx) (3.22a) 

a n p 

Q(x) = -f- (3.22b) 

a 2 

R(x) = — (3.22c) 

a 2 


Element matrices can then be derived using the method of weighted 
residuals in the same manner as previously described. In this case, 
element equations have the form, 


[[K c ] + [K h ]] < 


V 


= { Q c ) + { Q) 


(3.23) 


where the conductance matrices and heat load vectors are defined 
by 



(3.24a) 
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ty - 


Q {N} [Nj dx 


(3.24b) 


(Q c > 


- i 


- p ^ n 
dx 


(3.24c) 


{Q} = 1 R (N} dx 


(3 . 24d) 


Similarly, element equations for the two nodal unknowns are 


— — 




r -> 


r " 

*11 R 12 


T 1 




R io 


< 


” = (Q > + < 


> - T < 

? 

hi hi 


T 2 

c 

^2 

u 

*20 

_ - 




^ j 

1 

^ J 


(3.25) 


where K , i,j =0,1,2 is the summation of the corresponding 
coefficients in the conductance matrices [k ] and [k^] ; 


C* 2 dN. dN. r" 2 

hi - - p ^^ dx + 


x. 


Q N. N. dx 
13 


(3.26) 


i, j=0,l,2 


An additional advantage of using the self-adjoint differential 
equation is that the coefficients K^q and K^q shown on the right 
hand side of Eq. (3.25) are identically zero. This result can be 
proved by observing that the element interpolation function N^, 
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i = 1,2 are the solution of the homogeneous differential equation, 
Eq. (3.21), because is a linear combination of the homogeneous 

solutions f^(x) and f£(x) as s ^ lown Ec l s - (3.15a-b), i.e. 



dN t 

dx 


Q N. = 0 


i =1,2 


Multiplying this equation by the nodeless parameter interpolation 
function and performing integration by parts on the first term 

yields 


dN * 
_ 1 

dx 


X 2 *2 
r 


x. 


N, 


dN. dN 

P ~7~ ' dx + 

dx dx 


Q N. N. dx 

l 0 


= o 


Then since the nodeless interpolation function vanishes at 

nodes, i.e. at the coordinates x^ and x^, the above equations 
yield 


K i0 ■ 0 


i -1,2 


and the element equations, Eq. (3.25), become 


— 

— 






*11 

*12 


T i 

- - {Q c > 


\ 



< 


+ < 

> 

*12 

1 

CM 

CM 


T 2 

* J 



^2 


(3.27) 


After element nodal temperatures are computed, exact temperatures 
within an element can be obtained using the exact element temperature 
variation , Eq . (3.14). 
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To demonstrate the exact finite element formulation previously- 
derived, exact finite elements for eight heat transfer cases in 
several solids of different shapes and a flow passage (Fig. 4) are 
presented. In the first seven cases, heat transfer may consist of: 
(1) pure conduction, (2) conduction with internal heat generation, 

(3) conduction with surface heating, and (4) conduction with surface 
convection. Case eight is a one-dimensional flow where heat transfer 
may consist of fluid conduction and mass transport convection with 
surface heating or surface convection. For these cases, the boundary 
conditions considered are: 


or 


T = Constant 


(3.28a) 



(3.28b) 


or ‘ k ^ =h(T ‘ T «> ) (3.28c) 

where k is the material thermal conductivity, q is the specified 
surface heating rate per unit area, h is the convection coefficient, 
and is the convection medium temperature. In each case, the 

derivation of exact finite elements for appropriate heat transfer 
cases are given for clarity. Governing differential equations and 
the corresponding nodeless parameters, exact element interpolation 
functions, and element matrices for all cases are shown in Tables 1 
and 2 and Appendices A and E. 

3.2.1 Rod and Slab 


A rod element with arbitrary cross-sectional area A, circum- 
ferential perimeter p and length L as shown in (Fig. 4, Case 1) 



Case 1. ROD 


Case 2. SLAB 


V C 


A j 







V x -' / J 


\ y 


Case 3. HOLLOW CYLINDER Case 4. HOLLOW SPHERE 






Case 5. CYLINDRICAL SHELL Case 6. CONICAL SHELL 

\ , 

" L T 


flB 


Case 7, SPHERICAL SHELL Case 8, FLOW PASSAGE 


Fig. 4. Exact finite elements for one-dimensional conduction 
and convection cases. 




Table 1 


Governing Self-Adjoint Differential Equations 







Heat Loads 


Case 

Conduction 

(a) 

Convection 

(b) 

Convection 

(b) 

Source 

(c) 

Surface 

Flux 

(d) 


1 

--r-i 

dxLdx -1 

h£ T 
kA 

h£ T 

kA “ 

k 

32. 

kA 

2 

- — r — 1 

dxLdxJ 

— 

— 

2 

k 

— 

3 

d r dT-i 
dr L r dr-I 

— 

— 

k r 

— 

4 

dr 2dT"| 

drL r drJ 

— 

— 

b 2 

— 

5 

--r-i 

ds*-ds J 

-r^- T 
kt 

— T 
kt “ 

4 

k 

_9_ 

kt 

6 

--r^i 

dsL ds -1 

IT" sT 

kt 

h sT 

00 

kt 

k S 

s 

kt 

7 

dr s dTi 

dsL CO& a dsJ 

— 

— 

Q s 

i cos — 
k a 

q s 

cos 

kt a 

8 

_ _d[ P dll* 

dx<- dxJ 

JlE. P T 
kA 

7? P T 
kA 

— 

32. P 

kA 


*Combined conduction and mass transport convection where P = exp(-mcx/kA) 
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Table 2 

Nodeless Parameters for Thermal Problems 





T 0 


Case 

Convection 

(b) 

Source 

(c) 

Surface Flux 
(d) 


1 

T 


T 2 

ask. 

X 

oo 

2k 

2kA 

2 


aL 




2k 


3 






4kw 


4 


A. 




6k 



T 

QL 2 

qL 2 


JL 

oo 

2k 

2kt 

l : 

T 

i 

_Hbl 

O 

T co 

4kw 

4ktw 

7 

— 

5a! 

k 

5a 2 

kt 

8 

T 


qph 

oo 


me 


where w = ln(b/a) 



42 


is subjected to internal heat generation, surface heating, and 
surface convection. Governing differential equations for each heat 
transfer case in self-adjoint form are shown in Table 1. For example, 
the governing differential equation for the case of conduction with 
surface convection is 


" dl (kA + hpT = hpT “ (3.29) 

where k is the material thermal conductivity, h is the convection 
coefficient, and T^ is the convective medium temperature. A general 
solution to the above differential equation is 

T(x) = sinh mx + cosh mx + 

where m = ^hp/kA, and and C ^ are unknown constants. Applying 

the boundary conditions at the nodes, 

T (x=0) = T 1 and T(x=L) = T £ 

the two unknown constants are evaluated and the above solution 
becomes 


T(x) = (1 - 


sinh m(L-x) 
sinh mL 


sinh mx . 
sinh ml/ im 


. sinh m(L-x) . _ . sinh mx . _ 

^ sinh ml ; 1 'sinh ml/ 2 


(3.30) 


This exact element temperature variation can be written in the form 
of Eq. (3.14) where the element interpoation functions and the node- 
less parameter are: 
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N q (x) = (1 - 


sinh m(L-x) 
sinh mL 


sinh mx x _ 
sinh mL 9 0 


N x (x) 


sinh m(L-x) 
sinh mL * 


N 2 (x) 


sinh mx 
sinh mL 


(3.31) 


As described in the previous section, element equations for a 
typical self-adjoint differential equation have the form of Eq. (3.23), 
and using the definitions of the element matrices shown in Eq. (3.24), 
the element matrices for this problem are: 


L 

[k c ] =| kA {||} L-gj dx (3.32a) 

0 

L 

[ft h ] = 1 hp {N} [Nj dx (3.32b) 

L 

{Q} = [ hpT^ {N} dx (3.32c) 

J o 

where [K^] and [K^] are conductance matrices corresponding to 
conduction and convection, respectively, and {Q} is the load vector 
due to surface convection. With the exact interpolation functions 
shown in Eq. (3.31), the above element matrices can be evaluated in 
closed form. Exact nodal temperatures and element temperature varia- 
tion can then be computed using Eqs. (3.27) and (3.30), respectively. 

For the cases where the rod is subjected to an internal heat 
generation or specified surface heating, exact element interpolation 
functions and element matrices can be derived in the same manner as 
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described above. It should be noted that only the conductance matrix 
associated with conduction and heat load vectors corresponding to 
internal heat generation or surface heating exist in the two cases. 
The exact conductance matrix and heat load vectors are found to be 
identical to those from the conventional linear element. Therefore, 
exact nodal temperatures can also be obtained through the use of 
the conventional linear finite element in such cases. However, since 
the linear temperature variation is not an exact solution to these 
problems, the conventional linear finite element can not provide the 
exact temperature distribution within the element. 

The derivation of exact finite elements for one-dimensional heat 
transfer in a slab follows the derivation for the exact rod element. 

A slab with thickness L subjected to an internal heat generation 
(Fig. 4, Case 2) where both sides of slab may be subjected to a 
specified surface heating or surface convection. In Table 1, the 
governing differential equations are shown only for the case of pure 
conduction and conduction with internal heat generation because the 
effects of surface heating and surface convection enter the problem 
through the boundary conditions. For example, a governing differen- 
tial equation describing heat conduction in a slab with specified 
temperature T^ at x = 0 (node 1) and surface convection at x = L 
(node 2) is 


S ’- 0 


(3.33) 


where k denotes the material thermal conductivity. After solving 
for the general solution to the governing differential equation above 



and applying nodal temperatures as boundary conditions at x = 0 
and x - L, the exact element temperature variation is (see 
Appendix A) 
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T(x) 


(1 - f> T 1 + <f> T 2 - IV*> 



(3.34) 


With the corresponding element conductance matrix shown in Appendix B, 
exact element equations for this problem are 



j 


since at x = L (node 2) the boundary condition is 


dT (x=L) 
dx 


h(T 2 -T„) 


where h is the convection coefficient and is the surrounding 

medium temperature, therefore, the above element equations become 


r - 


r ^ 


r *v 

k k 


T 


. dT(x=0) 

L L 


1 


dx 


< 


» = * 


k k 

" L L +h 


T 2 


hT_ 

CO 


(3.35) 


L J l ^ ) 

the exact nodal unknown T 2 can then be computed and the exact 
element tempterature distribution is obtained using equation (3.34). 

The same procedure can be applied for the case when the slab is 
subjected to surface heating. In this case, the boundary condition 
is 
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q 


where q denotes the specified surface heating. When the slab 
consists several layers with different thermal conductivities, an 
exact element can be used to represent each layer. If the slab is 
subjected to surface heating or surface convection in addition, the 
above procedure applies for the elements located at the outer 
surfaces . 


3.2.2 Hollow Cylinder and Sphere 

A thermal model of a hollow cylinder with radial heat conduction 
subjected to an internal heat generation is shown in Fig. 4, Case 3. 
Specified heating or surface convection are considered through the 
boundary conditions at the inner and outer surfaces of radii a and 
b, respectively. Governing differential equations corresponding to 
each heat transfer case are provided in Table 1. For example, the 
governing differential equation for the case of pure conduction is 


k 


__d 

dr 


f 1 = ° 


(3.36) 


where k is the material thermal conductivity, and r is the radial 
coordinate. A general solution to the above differential equation is 

T(r) = C 1 + C 2 In r 


Nodal temperatures are imposed on the element boundary conditions. 


T(r =a) = T x 


T(r = b) = T 2 


and 
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and the exact element variation is obtained as (see Appendix A) , 

T ( r \ = | ln(b/r) ln(r/a) ■ 

^ ; *- w w J 

where w = ln(b/a). Note that the exact element variation for this 
case is completely different from the linear element variation, there- 
fore, the conventional linear finite element can not provide exact 
element or nodal temperatures. Applying the method of weighted 
residuals to the governing differential equation, element equations 
are 


r i 


(3.37) 


[K c ] (T } = {Q c } 


(3.38a) 


where 


(iy - f L f j r dr 

a 



< 

b 

= < 

. dT „ 
kr ~r— N 

> 


dr 




a 


(3.38b) 


(3.38c) 


Using the exact element interpolation functions shown in Eq. (3.37), 
element equations for this case are 



-ka 


kb 


dT 

dr 

dT 

dr 


(3.39) 
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When the cylinder is subjected to surface heating or surface 
convection, the same procedure previously described for the slab 
can be used. For example, in case of convection heat transfer on 
the outer surface, the boundary condition is 

r = b; " k S = h(T 2 - T »> 


where h is the convection coefficient and is the surrounding 

medium temperature. Thus, the element equations, Eq. (3.39), become 


k 

w 

_ k 
w 

< 

r 

H 

H 

v 

> = 

< 

r 'v 

. dT 

' ka dF 

> 

_ k 
w 

- + hb 
w 


T 

2 

>. j 



hb T 

oo 

L J 


(3.40) 


Exact finite elements can be formulated for conduction heat 
transfer in the radial direction of a hollow sphere with internal 
heat generation. A thermal model of a hollow sphere with inner and 
outer surface radii a and b, respectively, is illustrated in 
Fig. 4, Case 4. The hollow sphere may be subjected to surface 
heating or surface convection on both inner and outer surfaces. 

For heat conduction with internal heat generation, the governing 
differential equation is (see Table 1) 



t ' 2 §1 * <* 2 


(3.41) 


where k is the material thermal conductivity, Q is the heat 
generation rate per unit volume, and r is the independent variable 
representing the radial coordinate. A general solution to this 


differential equation is 
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^1 Qr^ 

T < r > =-T +c 2 -^k 


Due to the presence of the particular solution in the above general 
solution, a nodeless parameter exists, and the exact element 
variation is written in the form 


T (r) = N 0 (r) T Q + N^r) T ± + N 2 (r) T 2 (3.42a) 


where the element interpolation functions including the nodeless 
parameter are: 


N Q (r) = -^(r-a) (b-r) (r+a+b) ; T q = 6^ 


. , _ a(b-r) 

N l (r) " Tcb^y 


N 2 (r) “ r(b-a) 


(3.42b) 


Element matrices can be derived using the method of weighted residuals 
and element equations are resulted in the form 


[Kj {T} = {Q c } + {Q} 

where these element matrices are defined by: 


(3.43a) 




D 

■I 


, rdN-i idNi 2 , 

k {— } K— I r dr 
dr L dr J 


(3.43b) 


{Q c > - 


* 

s 

b 

, 2 dT .. 

k r -t— N 

l 

dr 



a 


(3.43c) 
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b 

{Q} = | Q{N} r 2 dr (3.43d) 

a 

If surface heating and surface convection are applied on the inner 
and outer surface, the same procedure described for the cylinder is 
required. 

3.2.3 Thin Shells 

Three thermal models of thin shells of revolution with cylin- 
drical, conical and spherical shapes are presented (see Fig. 4). 

These shells may be subjected to thermal loads such as surface 
heating, surface convection, and internal heat generation as shown 
in Fig. 4, Cases 5-7. In Case 5, a cylindrical shell of radius a, 
thickness t and meridional coordinate s is considered. Governing 
differential equations corresponding to different thermal loads are 
shown in Table 1. These governing differential equations are in the 
same form as for the rod element (Case 1). Therefore, the exact 
rod element interpolation functions and element matrices previously 
derived can be modified and used for the exact cylindrical shell 
element . 

A truncated conical shell element with thickness t is shown 
in Fig. 4, Case 6. Governing differential equations corresponding to 
internal heat generation and surface heating are given in Table 1. 
These differential equations are in the same form as for the hollow 
cylinder (Case 3) with surface heating, and therefore, element 
interpolation functions and element matrices are similars. For 
the case of the shell subjected to surface convection, a form of 
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nonhomogeneous modified Bessel's differential equation results. 


d 2 T + ldT__h T = __h_ T 
, 2 s ds kt kt 00 


(3.44) 


A general solution to the above differential equation includes 
modified Bessel functions of the first and second kind of order zero. 
A nodeless parameter also exists in this case due to the nonhomo- 
geneous differential equation. Applying nodal temperatures as the 
boundary conditions at s = a and s - b, exact element interpola- 
tion functions are obtained as shown in Appendix A. 

Fig. 4, Case 7 shows a truncated spherical shell with radius 
a and thickness t. The spherical shell may be subjected to 
internal heat generation or surface heating. Governing differential 
equations corresponding to these thermal loads are in the form of 
Legendre f s differential equation of order zero. For example, the 
governing differential equation for the case of uniform surface 
heating q is 


(1 - 


n 2 ) 


d 2 T 

dn 2 


. 2 „ - - ¥r 


dn 


kt 


(3.45) 


where q = sin (s/a). A general solution to the above differential 
equation is 


T 


2 


In [ 


1+n 

1-n 


J + c 


2 



(3.46) 


where and C ^ are unknown constants. By imposing nodal 

temperatures as element boundary conditions at s = 0 and s = L, 
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exact element interpolation functions are obtained as shown in 
Appendix A. 

Due to the complexity of the exact element interpolation func- 
tions that arise from the truncated conical shell with surface 
convection and the truncated spherical shell, the corresponding 
element matrices in closed form are not provided. The element 
matrices, if desired, can be obtained using the element matrix 
formulation shown in Equations (3.14a-d) and performing the integra- 
tions numerically. 

3.2.4 Flow Passage 

A thermal model of fluid flow in a passage with conduction and 
mass transport convection is illustrated in Fig. 4, Case 8. The 
fluid may be heated by surface heating, or surface convection. 
Governing differential equations corresponding to these heat transfer 
cases are given in Table 1. For simplicity, consider the case without 
heat loads where the governing homogeneous differential equation is 
given by 


1 [kA —1 

dx L dx 


. • d T _ 
+ me — = 0 
dx 


(3.47) 


where k is the fluid thermal conductivity, A is the flow cross- 
sectional area, m is the fluid mass flow rate, and c is the fluid 
specific heat. A general solution to this differential equation is 

T(x) = c i + C 2 ex P (2“x) 

where and are arbitrary constants and a = mc/2kA. An 

exact finite element with length L and nodal temperatures and 
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at x * 0 and x « L, respectively, can be formulated. The 
exact element temperature variation is 


T(x) = [l “ “ 
1 


2ax 
e 

2aL 


e 



(3.48a) 


As previously described, the appearance of the first-order derivative 
term in the governing differential equation results in an unsymmetrical 
conductance matrix (see Eq. (3.18b)). In this case, the corresponding 
element conductance matrices are 


[Kj 


kAa 


, 2otL 
(e 

, 2aL 


+ 1 ) 

- 1 ) 


1 

-1 



' K vl 


me 

2 


-1 

-1 


1 

1 


(3.48b) 


(3.48c) 


where [K ] and [K ] denote conductance matrices representing 
c v 

fluid conduction and mass transport fluid convection, respectively. 

It has been shown that if the conventional finite element with 
an optinum upwind weighting function is used, exact temperatures at 
nodes can also be obtained [21 3 . With upwind weighting functions 
the element temperature variation is expressed as. 


T (x) = |l - f + F(x) 



(3.49a) 


where F(x) is the optimum upwind weighting function defined by 
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F(x) = [coth (2aL) - [3(^ - f )] 

L 


With these element interpolation functions, element conductance 
matrices corresponding to the fluid conduction and mass transport 
convection are 


[K ] 

c upwind 


kA 

L 


1 -1 


-1 1 


(3.49b) 


[K*] 


upwind 



‘-1 l" 



~ l -l" 

me 

2 


• 

+ — (coth (aL) 

9 




1 

M 

H* 

L_ 

Z 


1 

H* 


(3.49c) 


It can be shown that the combination of these element conductance 
matrices are identical to those obtained from the exact finite 
element, Eqs. (3.48b-c). Therefore, the conventional finite element 
with the optimum upwind weighting function provide exact nodal 
temperatures. However, since the upwind element temperature varia- 
tion differs from the exact element temperature variation shown in 
Eq. (3.48a), the finite element with the optimum upwind weighting 
function does not provide the exact temperature variation within an 
element. 


3.3 Exact Finite Elements in Thermal-Structural Problems 

With the general exact finite element formulation described 
in section 3.1, exact structural finite elements can be developed 
for problems governed by ordinary differential equations. For 
example, exact finite elements for a rod loaded by its own weight 


55 


or a beam with a distributed load can be formulated. However, for 
the purpose of demonstrating benefits on exact finite elements in 
coupled thermal-structural problems, exact structural finite elements 
subjected to thermal loads are considered herein. 

3.3.1 Truss 

Typical thermal and structural models for truss elements are 
shown in Fig. 5. For a steady-state analysis, exact thermal finite 
elements for internal heat generation, surface convection and 
specified surface heating are presented in section 3.2. In this 
section the exact element temperatures are used in the development 
of truss elements for computations of displacements and thermal 
stresses. 

For a truss element subjected to a temperature change, thermal 
strain is introduced in the stress-strain relation; 

a x = E [ ^r~ " a(T(x) “ T ref>] < 3 - 5 °) 

where a x is the axial stress, E is the modulus of elasticity, 
u is the axial displacement which varies with the axial coordinate 
x, a is the coefficient of thermal expansion, T(x) is the 
temperature, and T re f is the reference temperature for zero stress. 
The rod equilibrium equation with an assumption of negligible body 
force is 



(3.51) 


which when combined with the stress-strain relation, Eq. (3.50), and 
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Fig. 5. Thermal and stress models of rod element. 
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multiplied through by the truss cross-sectional area A, yields 
the governing differential equation. 


EA = aEA ^ (3.52) 

dx 2 dx 

Since the temperature T is known from the thermal analysis, a 
general solution to the above differential equation can be obtained • 

An exact finite element can be formulated by applying the nodal 
displacements u^ and u 2 as the boundary conditions at x = 0 
and x = L, respectively. In this case, the exact element displace- 
ment variation is 


x L 

u(x) - (o J T dx - a £ j* T dx) + (1 - •£) U]L + (£) u 2 (3.53) 

0 0 


or in the matrix form 


u (x) = [NqCx) N 1 (x) N 2 (x)J < 




u„ 


u. 


- w (“> 


(3.54a) 


where N^(x) is the element nodeless interpolation function; 

N^, i “1,2 are typical element interpolation functions, u^ is 
the nodeless parameter, and u^, i =1,2 are the element nodal 
displacements. The element interpolation functions are 


x L 

N q (x) = a j T dx - a ^ j 
0 0 


T dx 


(3.54b) 
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N 1 (x) = 1 - f N 2 (x) = f (3.54c) 

where, for convenience, the nodeless parameter is taken as unity 

in this case. Note that the element nodeless interpolation function, 
Nq(x), vanishes at nodes and depends on the integrals of element 
temperature variation obtained from the thermal analysis. 

To derive exact element matrices, the method of weighted resid- 
uals is applied to the equilibrium equation (3.51). After performing 
an integration by parts and using the stress-strain relation, 

Eq. (3.50), element equations and element matrices are obtained. 

These element equations are in the same form as those obtained from 
the variational principle described in section 2.3 and can be 
expressed as 


[KJ {u} = {F x } (3.55) 

where [K ] is the structural element stiffness matrix, {u} is 
s 

the vector of nodal displacements, and {F^} is the equivalent 
nodal thermal load vector. The element matrices are defined by (see 
Eqs. (2.27a) and (2.27d)) 

r L dN dN 

[ V - A*J L-jfJ “X < 3 - 56 “> 

0 

f L dN 

= AEa ) { lf } (T - T ref> dx 
0 


(3.56b) 
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Using the exact displacement interpolation functions, Eq. (3.54a), 
the element stiffness matrix above is a three by three matrix which 
contains coefficients K_ , i,j =0,1,2. Since the governing 
differential equation, Eq. (3.52), can be cast in the self-adjoint 
form (see section 3.2), this element stiffness matrix is symmetric 
and i =1,2 are zero. Both the element stiffness matrix and 

the equivalent nodal thermal load vector can be evaluated in closed 
form as, 




i -l 


-l 




-l 




(3.57a) 


(3.57b) 


where 


L 

F t = AEa | (T - T ref ) dx (3.57c) 

0 

Once exact nodal displacements are determined, exact displacement 
variation within an element can be computed from Eq. (3.53). Exact 
element stress can also be obtained by substituting element displace- 
ment variation, Eq. (3.53), into the stress-strain relation, 

Eq. (3.50). In this case, the exact element stress in terms of 


nodal displacements is 
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a 

x 



ot 

L 


r 


< T - T ref> 



(3.58) 


Using the exact element temperature variations obtained from the 
thermal analysis (cases la-ld) , both element nodeless interpolation 
functions N^(x), Eq. (3.54b), and the equivalent nodal thermal load 
F , Eq. (3.57c), can be evaluated in closed form as shown in Table 
3 and Appendix B, respectively. 


3 .‘3 . 2 Hollow Cylinder 

For a hollow cylinder where the temperature T varies only in 
the radial direction (Fig. 6), the only non-zero displacement is 
u(r) and all shearing stresses are zero. The radial stress a 
and circumferential stress cf^ satisfy the equilibrium equation [22] 


da 


r 


dr 



= 0 


(3.59) 


The stress-strain relations are 

£ r [<7 r - v(a 9 + + a < T - T ref> 

% = E [0 e " V(0 r + °z )] +a(T ‘ T ref ) 

£ z = E [a z " v (a r + a 0^ + a (T " T ref ) 

where v is Poisson's ratio; e , £q and e z are 


(3.60a) 


(3.60b) 


(3.60c) 


the radial. 


circumferential and longitudinal strain, respectively. 
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Table 3 

Truss Element Displacement Interpolation Functions, N S (X)* 


Case 


N 0 (X) 


1(a) 


Ot(T« T-.)L ry &T-.L ry « 

— - — (X -X) + — — (-X + 3X - 2X ; 

2 6 


1(b) 


T -T cosh mL + T n (cosh mL-1) 

a ( [(cosh mLX -1) 

m sinh mL 


T -T 

- X(cosh mL-1)] *4 - (sinh mLX - X sinh mL) } 


m 


1(c) 


“< I 2- T l )L ,„2 . “V 


(X -X) + — — (-X + 3X - 


2X 3 ) 


1(d) 


2 aT 0 L 2 
=-r -± — (X -X) + — (-X + 3X 


2X 3 ) 


*For all cases: N^X) = 1-X, N 2 (X) = X where X « x/L. 
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ONE DIMENSIONAL THERMAL MODEL 

Fig. 6. Thermal and stress models of axisymmetric element. 
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For the case of a thin hollow cylinder, the assumption of plane 

stress (a_ = 0) is used. Substituting the stress-strain relations, 
z 

Eqs. (3.60a-b), into the equilibrium Eq. (3.59) and using the strain- 
displacement relations. 

e = 4^ and = — (3.61) 

r dr or 

where u denotes the radial displacement, the governing differential 
equation for the case of plane stress is 


d r_l d(ru) -I 
dr L r dr 1 


(1 + 


. dT 

v) “-37 


(3.62) 


A general solution to this* differential equation is given by 

r 

u(r) * (1 + v) y j (T - T ref ) r dr + C x r + (3.63) 

0 


Since the radial temperature variation T is known from the 
thermal analysis (see section 3.2.2), the exact axi symmetric element 
displacement variation can be derived by applying the nodal displace- 
ments u^ and as the boundary conditions at r = a and r = b, 

respectively. The exact element displacement variation is 


u(r) - (1 + v) ^ (T - T ref ) v dr 

b 

S 


- a + .)|4H 

r (b^ - a 1 


(T - T ref ) r dr 
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+ 


,, 2 2 N 
r a (b ~ r ) •) 

r (b 2 - a 2 ) J 


U 1 + 



(r 2 - 
(b 2 - 



(3.64) 


or in the matrix form 


u(r) = [NgCir) N 1 (r) N 2 (r)J 


« u. 


“ lAJ 


(3.65a) 


where N^Cr) is t ^ ie element nodeless interpolation function; ISL , 
i =1,2, are typical element interpolation functions, u^ is the 
nodeless parameter, and u^, i =1,2 are the element nodal displace- 
ments. The element interpolation functions are 


r 

N Q (r) = (1 + v) 2L j (x - T ref ) r dr 

a 


- (1 + v) 


a (r 2 - a 2 ) 
r (b 2 - a 2 ) 


b 



a 


(3.65b) 


N x (r) 


? ? 

r a. (b Z ~ r“) , 
r (b 2 - a 2 ) 


and 


N 2 (r) 



(r 2 - 


a 2 ) 


(b 2 - 


* 2 ) 


(3.65c) 


Like for the exact truss element, the nodeless parameter u^ is 
taken as unity, and the element nodeless parameter Ng(r) vanishes 
at nodes. Element matrices can be derived by following the same 
procedure described for the truss element. In this case, the element 
stiffness matrix and the equivalent nodal thermal load vector are 
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b 



[B S ] T [d] 


[B s ] r dr 


(3.66a) 


b 

{F T } = J [B s ]T {a} (T " T ref ) r dr (3.66b) 

a 


where [B g ] is the strain-displacement matrix obtained from 
Eq. (3.61), 




du 


“dN ] _ 

dN2 ~ 
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£ 

r 

< 

t = < 

dr 

► = 

dr 

dr 

< 

-L 

► = [B s ] {u} 

c 


_u 


Nl 

^2 


U 2 


£ e 

k 4 


r 


r 

r 


^ - 



[D] is the elasticity matrix (plane stress), 


(3.67a) 


[D] 




1 


(3.67b) 


and (a) is the vector of coefficients of thermal expansion, 


1 


{a} = 


a 1 > 


1 


(3.67c) 


Using the exact element interpolation functions, Eq. (3.65), 
the element stiffness matrix and equivalent nodal load vectors 
corresponding to the heat transfer cases (see section 3.2.2) can be 
derived in closed form. Due to complexity of the element interpolation 
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functions, a computer-based symbolic manipulation language MACSYMA 
was used to perform the algebra and calculus required for these 
element matrix derivations. Results of these element matrices and 
exact element interpolation functions are shown in Appendix B and 
Table 4, respectively. 

Once nodal displacements are computed, exact thermal stresses 
in both radial and circumferential directions can be determined. 
Using the stress-strain relations, Eq. (3.60), and the strain- 
displacement equations in the form of Eq. (3.67a), the element 
stresses can be written in terms of nodal displacements as 


* y 

o 

r 

r 


< 

* - [D] « 

°Q 

. ' 4 



[B s ] {u} - (1 + v) o(T - T ref 


,r 


(3.68) 


For the plane strain case (c z = 0), all equations formulated 

for the case of plane stress above may be used by replacing 
2 

E.(l-v ) for E, v/(l-v) for v, and (l+v)a for a. In addi- 
tion, the longitudinal stress exists in this case and can be computed 
from the last equation of the stress-strain relations, Eq. (3.60c), 

a z = v(a r + a 0 } ~ aE (T “ T ref } (3.69) 


3.4 Applications 

To demonstrate the capabilities of the exact thermal and 
structural finite elements developed in sections 3.2 -3.3, the finite 
element thermal analysis program TAP2 [23] and the finite structual 
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Table 4 

Axisyrametric Element Displacement Interpolation Functions 


2 2 2 2 
Z 2. ,bv (b -r )a w 


V r) ■ 4^ (<T 1 + 72 "V I r ^ - 
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,,2 2 . 

(b -a ) 


+ (i, + wt 0 > [r 2 i„(f) - (r2 -; 2 >f» : 

Z 0 a (b -a Z ) 


„ . 2 2 W 2 2, 

+ v 2 ; (b -c Hr -a ) )} 

2b 


2 2 

r rn a(b Z -r Z ) 

V r) ' Z2 2. 

r (b -a ) 


2 2, 

N,<r) = b(r 
1 r(b -a Z ) 
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analysis program STAP [24] are used. Elements discussed in this 
chapter were added to these programs. Conventional finite elements are 
also available in these programs, so comparisons between exact finite 
elements and conventional finite elements could be made. 

3.4.1 Coffee Spoon Problem 

The exact rod element for conduction and convection described 
in section 3.2.1 is used for one-dimensional heat transfer in a 
coffee spoon [25], Fig. 7. The lower-half of the spoon submerged in 
coffee is convectively heated by the coffee at 339 K, and* the upper- 
half is convectively cooled by the atmosphere at a temperature of 
283 K . The ends of the one-dimensional spoon model are assumed to 
have negligible heat transfer. 

Three finite element models are used to represent the spoon: 

(1) two exact finite elements, (2) two linear conventional finite 
elements, and (3) ten linear conventional finite elements. Tempera- 
ture variations computed by these three finite element models are 
compared in Fig. 7. The figure shows that two conventional finite 
elements predict nodal temperatures with fair accuracy but are unable 
to provide details of the nonuniform temperature distribution includ- 
ing the zero temperature gradients at both ends of the spoon. The 
temperature variation obtained from ten conventional finite elements 
is in excellent agreement with the result from two exact finite 
elements. It should be noted, however, that an approximate solution 
results from the use of conventional finite elements since the exact 
solution to the problem is in terms of hyperbolic functions, which 
were used in the exact element interpolation function (Appendix A, 


Case lb) . 
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Fig. 7. Conventional and exact finite element solutions 
for coffee spoon with conduction and convection. 
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3.4.2 Thermal Stresses in Hypersonic Wing 

A 136 member truss model of a hypersonic wing [26], Fig. 3, was 
chosen to illustrate the use of exact truss finite elements. The wing 
is assumed to have varying convective heat along the leading edge, top 
and bottom surfaces and is convectively cooled internally. Tempera- 
tures along the wing root are specified. Two finite element thermal 
models are used to represent the wing truss. The first model 
consists of 136 exact condution-convection rod elements (see 
section 3.2.1) with one element per truss member. The second model 
is identical to the first model, but linear conventional finite 
elements are used. Fig. 9 shows a comparison of temperature distri- 
butions along the bottom members of the center rib of the wing truss. 
Results show that the exact finite element model provides a realistic 
temperature distribution which is characterized by higher temperatures 
near the center of each truss member and lower temperatures at the 
nodes. The conventional finite element model underestimates the 
actual temperatures and is not capable of capturing the highly 
nonlinear temperature distribution along the rib. Therefore, further 
mesh refinement of the conventional finite element model is needed 
if a realistic temperature distribution is to be predicted. 

For the structural analysis, both models employ the same 
discretization as in the thermal analysis. The structural boundary 
conditions consist of constraining the nodes along the wind root. 

Truss member temperatures obtained from the exact finite element 
thermal model are directly transferred to the exact finite element 
structural model, for computations of displacements and thermal 
stresses (see section 3.3.1). Likewise, displacements and thermal 




Fig. 8. Thermal structural truss model of a hypersonic wing. 
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stresses computed from the conventional finite element structural 
model are based upon linear member temperatures obtained from the 
conventional finite element thermal model. Comparison of the thermal 
stress distributions for the two analyses are made as shown in 
Fig. 9. The figure shows that conventional finite elements under- 
estimate member stresses with a relatively large error. This error 
is caused by the use of the inaccurate temperature distribution 
from the conventional finite element thermal model. Comparative 
temperature and stress distributions of other wing sections (not 
shown) have similar trends. The results clearly demonstrate that 
improved thermal-structural solutions can be obtained through the 
use of exact finite elements. 


Chapter 4 


MODIFICATION OF EXACT FINITE ELEMENT FORMULATION FOR 
ONE-DIMENSIONAL LINEAR TRANSIENT PROBLEMS 

In the preceding chapter, exact thermal finite elements for one- 
dimensional steady-state heat transfer problems were presented. 
Steady-state element temperature interpolation functions were 
formulated in closed form based upon solving ordinary differential 
equations. In transient analysis, exact element temperature inter- 
polation functions cannot be obtained in closed form since general 
solutions to typical transient problems are infinite series. However, 
by modifying the steady-state element temperature interpolation 
functions for the transient analysis, improved transient temperature 
solutions can be obtained as described in this chapter. 

In steady-state analysis, finite element temperature distribu- 
tions are a function of only the spatial coordinate, but for transient 
analysis, the element temperature distribution are a function of both 
space and time. For example, a one-dimensional transient heat conduc- 
tion is governed by the partial differential equation. 


, . 3 2 T . 3T 

kA j = pcA — 

3x 9t 


(4.1) 


where k is the material thermal conductivity, p is the density, 
c is the specific heat, A is the conduction area, and T is the 
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temperature which varies with the spatial coordinate x and time t. 
A two-node linear conventional finite element may be used in the 
analysis where the element temperature variation T is expressed 
in the form (Fig. 10(a)) 


T(x,t) - |1 - L lJ i 


T x (t) 


T 2 (t) 


> = [^(x) N 2 (x)J < 


T L (t) 


T 2 (t) 


(4.2) 


where N^(x), i — 1,2 are the element interpolation functions which 
are a function of the spatial coordinate x; L is the element length, 
and T\(t), i =1,2 are the time-dependent nodal temperatures. 

With the heat equation shown in Eq. (4.1), the corresponding 
element equations and element matrices can be derived as described 
in section 2.3. Typical element equations have the form 


[C] {T} + [K] {T} = {Q} (4.3) 

where {T} and (f) denote vectors of nodal temperatures and the 
time rate of change of nodal temperatures, respectively. The matrix 
[K] and the vector {Q} represent the conductance matrix and the heat 
load vector, respectively, and have the same meaning as previously 
described for the steady-state analysis in the preceding chapter. 

The additional matrix [c] is called the capacitance matrix and defined 
by (see Eq. (2.15a)) 


L 

[C] = j* pcA {N t > [N t ] dx 
0 


(4.4) 



(a) Conventional F. E. 


x 




Fig. 10. One-dimensional element interpolation functions 



For the linear interpolation functions shown in Eq. (4.2), the 
capacitance matrix can be evaluated as 
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[C] 



1 

2 


(4.5) 


This form of the capacitance matrix, Eq. (4.5), is called a consistent 
capacitance matrix because its definition is consistent with the 
matrix formulation, Eq. (4.4). Quite often, the above capacitance 
matrix is approximated by lumping the off-diagonal terms with the 
diagonal terras to give, 


[c] = 


l o 

0 1 


(4.6) 


and is called a lumped capacitance matrix. It should be noted that 
degradation of the solution accuracy may result from the use of the 
lumped capacitance matrix compared with the consistent capacitance 
matrix. However, computational advantages (e.g. explicit time 
integration algorithms) may be achieved using the lumped capacitance 
matrix whereas the loss of solution accuracy may be insignificant 

[5]. 


4.1 The Nodeless Variable 

In the preceding chapter, exact finite elements for steady- 
state analysis are formulated based upon solving ordinary differential 
equations. Exact element temperature variations after imposing nodal 
temperatures as boundary conditions are written in the form, 
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T (x) = N q (x) T q + N x (x) T 1 + N 2 (x) T 2 (4.7) 

where N^, i =1,2 are element interpolation functions; T^, i =1,2 
are unknown nodal temperatures, Nq(x) is the element nodeless 
interpolation function, and Tq is a known nodeless parameter. 

For transient thermal problems, it is not possible to formulate 
exact element interpolation functions in closed form because general 
solutions to typical transient problems are infinite series. However, 
since the transient response may approach exact steady-state solutions 
as time becomes large, the use of the exact steady-state element 
temperature variation in the form of Eq. (4.7) may provide better 
accuracy of solutions than those obtained from the linear conventional 
element, Eq. (4.2). 

To use the steady-state element temperature variation for 
transient analysis, Eq. (4.7) is written in the form, 

T(x,t) = Nq(x) T q + N x (x) T x (t) + N 2 (x) T 2 (t) (4.8) 

where the unknown nodal temperatures T^ and T 2 become a function 
of time t. Since the nodeless parameter Tq is known and independent 
of time, the product of the nodeless interpolation function and the 
nodeless parameter, Nq(x) Tq, retains the same shape throughout 
the transient response. Characteristics of the element temperature 
variation expressed by Eq. (4.8) during the response are illustrated 
in Fig. 10(b) . 

Equation (4.8) may not be a good representation for a transient 
thermal response as will be shown by the following argument. As 
described in section 3.1.1, the nodeless parameter Tq 


is a scalar 
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quantity which contains a physical parameter associated with a given 
heat load- For example, the nodeless parameter for a slab subjected 
to a uniform internal heat generation rate is given by (see 

Table 2) 



where k is the material thermal conductivity and L is the 
element length- If a slab is modeled by an exact finite element, 
the element temperature variation is given by (see Appendix A, 
Case 2) 


- (i - 
L V I / 


QjL 

2k 


a -f) 


T i + 


(— ) T 

V 2 


where and are the element nodal temperatures. If both 

surfaces of the slab have a specified temperature T g in addition, 
the above equation becomes 


T(x,t-0) =f (1 ~ f) 



+ - ?> T s + ( f> T s 


(4.9) 


For the case where the internal heat generation is raised instan- 
taneously from to Q^, the transient temperature variation 

within the slab should gradually increase and reaches the new steady- 
state temperature variation 


Q 0 L 

T(x,t-~) “ f (1 - f) 2 


21c 


+ (1 - T + (f) T 

L s L S 


(4.10) 
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where the new nodeless parameter is 



Since the nodal temperatures at both sides of the slab are fixed, 
the above two temperature variations, Eqs. (4.9 -4.10), suggest that 
the nodeless parameter Tq should vary with time so the element 
temperature can change gradually during the transient response. This 
argument leads to a modification of the element temperature variation 
employed in the steady-state analysis, Eq. (4.8), to the form 

T(x,t) = N q (x) T Q (t) + (x) T x (t) + N 2 (x) T 2 (t) (4.11) 

where the nodeless parameter becomes an additional time-dependent 
element unknown and is called a nodeless variable. A typical element 
temperature variation with the nodeless variable Tg(t) is illus- 
trated in Fig. 10(c). 


4.2 Element Equations and Matrices 


In this section, element equations and matrices for both the 
nodeless parameter approach and nodeless variable approach are 
presented. In the nodeless parameter approach, the element tempera- 
ture variation shown in Eq. (4.8) can be written in the matrix form 



T x (t) 

T 2 (t) 


T(x,t) = LN 0 (x) N^x) N 2 (x)J 


< 


> 


(4.12) 
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Using the definition of the element capacitance matrix shown in 
Eq. (4*4), the capacitance matrix for the given element interpolation 
functions can be derived. Element equations, Eq. (4.3), can then be 
written explicitly as 


o 

o 

o 

C 01 

C 02 


3T Q /3t 


00 

0 

0 

C 01 

C 11 

C 12 

-< 

3^/at 

► + 

0 

K 11 

K 12 

C 02 

C 12 

C 22 


3T 2 /3t 


0 

K 12 

K 22 


(4.13) 


where c -j_j » =0,1,2 are typical terms in the capacitance matrix; 

and Q^, i,j =0,1,2 are typical terms in the element stiffness 
matrix and the heat load vector previously described in the steady- 
state analysis. Since the nodeless parameter Tq is constant, its 
time-derivative 9 Tq/ 3t is zero. Therefore, the first equation 
which involves the nodeless parameter is uncoupled from the nodal 
unknowns in the second and third equations. Hence, the above equa- 
tions reduce to 
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(4.14) 


Note that these equations contain two basic element nodal unknowns 
as for the linear conventional finite element. Once the nodal 
temperatures at a typical time are computed, element temperature 
variation can be obtained using Eq. (4.12). 

In the nodeless variable approach, the element temperature varia- 
tion shown in Eq. (4.11) can be written in the matrix form 
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T(x,t) = L n q( x ) ^(x) 


n 2 (x)J 


T 0 (t) 

< T L (t) > 


T 2 (t) 


( 4 . 15 ) 


where Tq( 0 denotes the element nodeless variable which is a func- 
tion of time as the unknown nodal temperatures T^(t) and T^Ct), 
Since the element interpolation functions are identical to those 
used in the nodeless parameter approach, the element matrices are 
also identical. Element equations obtained using this approach have 
the form 


TP 

o 

o 

C 01 

C 02 


' 

*0 

o 

o 

h~* 

C 11 

C 12 

< 

*1 > 

C 02 

C 12 

C 22_ 


T 

i 2 



( 4 . 16 ) 


Because the nodeless variable is unknown, the equations are coupled 
through the capacitance matrix due to the presence of Tq. Thus 
typical element equations obtained from the nodeless variable approach 
contain three unknowns, i.e. one more unknown than the nodeless 
parameter apporach or the linear conventional finite element. 

An advantage of the nodeless parameter and nodeless variable 
approaches is that both can provide an exact steady-state solution 
at the initial condition for the transient response. As time becomes 
large and new steady-state thermal equilibrium is reached, the exact 
temperature distribution may be predicted by both approaches. It 
should be noted that with the use of the nodeless variable approach 
the temperature variation within the nodeless variable element can 
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vary with time even though the element nodal temperatures are fixed. 
This feature is characterized by the term N^(x) TqCO shown in 
Eq. (4.11) and is different from the linear conventional finite 
element where the element temperature distribution is completely 
controlled by nodal temperatures. 


4.2.1 Rod Element 

The rod element with heat conduction combined with surface 
convection, internal heat generation, or surface heating previously 
considered in Fig. 4, Cases la-d is extended for transient analysis. 
For each heat transfer case, the governing differential equation for 
the temperature distribution T(x,t) can be derived using an energy 
balance on a small segment of the rod. These governing differential 
equations are: 

pcA - kA =0 (4.17a) 

3t 3x 2 


pcA - kA + hpT = hpT^ (4.17b) 

3t 3x Z 


PcA -r~ - kA = QA 

3t 3x 2 


(4.17c) 


. 3T . . 
pcA — - kA 


3 2 T 

3x 2 


= qp 


(4. 17d) 


where A is the element cross-section area, h is the convection 
coefficient, p is the cross-section perimeter, T m is the surround- 
ing medium temperature, and q is the specified surface heating rate 
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per unit area. As one example, conduction with internal heat genera- 
tion where the exact steady-state element temperature variation is 
(see Appendix A, Case 1(b)) 


TCx.t) - z (1 - £> T 0 + (1 - f) Tj + (f> T 2 


f 


- [N q (x) N 3_(x) N 2 (x)J \ 


T 1 > 


(4.18a) 


where Tq is the nodeless parameter given by (see Table 2) 


T 

0 2k 


(4.18b) 


Using the exact element interpolation functions shown in Eq. (4.18a) 
above, the capacitance matrix is derived using Eq. (4.4.). The 
conductance matrix and heat load vector are derived using Eqs. (2.15b) 
and (2.16b), respectively. Therefore, the element equations are 
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(4.19) 


In the nodeless parameter approach, the nodeless parameter Tq 


is constant and the above element equations reduce to 
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with two unknown nodal temperatures T^(t) and T^Ct)* It should 
be noted that the element equations obtained from using the nodeless 
parameter approach shown in Eq. (4.20) above are identical to those 
obtained from the linear conventional finite element for this heat 
transfer case. Thus, results of nodal temperatures during the 
transient response are also identical. However, results of element 
temperatures are different due to the difference of their element 
interpolation functions, Eqs. (4.2) and (4.18b). As the transient 
response reaches the steady-state, the nodeless parameter approach 
provides exact solution for both nodal temperatures and element 
temperature variations where only exact nodal temperatures are 
obtained through the use of the linear conventional finite element. 

In the nodeless variable apporach, the element equations with 
two unknowns of nodal temperatures and an unknown of nodeless 
variable shown in Eq. (4.19) must be solved simultaneously. It can 
be seen from these equations that as the steady-state thermal 
equilibrium is reached, the rate of change of nodal temperatures 
and the nodeless variable vanish. Then the first equation yields 


T 


0 



which is identical to the nodeless parameter shown in Eq. (4.18b). 
This means the nodeless variable varies during the transient response 
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and provides the value required for computation of the exact 
temperature variation when thermal equilibrium is reached. 

For other heat transfer cases such as conduction with surface 
convection or surface heating, the same procedure is applied. 

Element matrices corresponding . to each heat transfer case in the form 
of Eq. (4.16) are given in Appendix C. Capabilities of the nodeless 
parameter and the nodeless variable finite elements for transient 
analysis are evaluated by comparisons with an exact transient 
conduction-convection solution and the linear conventional finite 
element in the first example at the end of the chapter. 


4.2.2 Axisymmetric Element 

Similar to the rod element, the axisymmetric element previously 
described in the steady-state heat transfer (Fig. 4, Case 3) is 
extended for the transient analysis. Radial heat conduction is 
combined with internal heat generation and specified surface heating 
or surface convection on the inner or outer cylinder surfaces are 
considered through the boundary conditions. The governing differ- 
ential equations for the cases of pure conduction and conduction 
combined with internal heat generation are 


pc 


3T 

at 


k JL 

r 9r 
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(4.21a) 
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3T 
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k a , 3Tv - 

7 77 (r 77> * Q 


(4.21b) 


respectively, where r denotes the radial coordinate. 

Element equations can be derived by the method of weighted 
residuals applied to the governing differential equations. 
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Typical element equations are in the form of equation (4.16). The 
element conductance matrix and heat load vector are identical to those 
obtained in the steady-state analysis shown in Eqs. (3.38b) and 
(3.38c), respectively. The element capacitance matrix associated with 
the rate of change of nodal temperatures has the form 

b 

[C] = j pc {N t > [N t J r dr (4.22) 

a 

For the element interpolation functions associated with the heat 
transfer cases shown in Appendix A, the corresponding capacitance 
matrix can be evaluated in closed form. Capacitance matrices in 
the form of Eq. (4.16) are given in Appendix C. 

4.3 Applications 

4.3.1 Transient Heat Conduction in a Rod with Surface Convection 
A rod with length L subjected to surface convection and 
specified end temperatures is shown in Fig. 11(a). Initially the 
rod is convectively cooled by a surrounding temperature at 255 K 
and, at time t = 0 , the surrounding temperature is raised 
instantaneously to 589 K. Transient temperature distributions along 
the rod are computed using: (1) the exact solution [27], (2) two 

linear conventional finite elements, (3) two nodeless parameter 
finite elements, and (4) two nodeless variable finite elements. In 
each finite element model, the element lengths are taken to be equal 
(L/2) with an unknown of nodal temperature at the center of the rod. 
Comparisons of the temperature variations obtained from these three 
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finite element models and the exact solution for t = 0, 0.01 and 

0.3 s. are made as shown in Fig. ll(b-d). 

At time t * 0 while the rod is in thermal equilibrium, two 
nodeless variable finite elements provide the exact steady-state 
temperature distribution. Two conventional finite elements are 
unable to provide details of the nonuniform temperature distribution 
due to the assumption of linear temperature distribution within the 
element. At time t * 0.01 s. (Fig. 11(c)) after the rod has been 
convectively heated, the differences in the transient response 
predicted by three finite element models are shown clearly. Two 
linear conventional finite elements predict the unknown nodal 
temperature at the center of the rod with fair accuracy but the 
element temperature distributions are overestimated from the actual 
temperature distribution with a relatively high error. Two nodeless 
parameter finite elements yield the unknown nodal temperature with 
the same accuracy as of two linear conventional finite elements but 
predict extremely poor element temperature distributions. Two nodeless 
variable finite elements provide the best approximation of the 
unknown nodal temperature with excellent temperature distributions 
within the elements. As the rod temperatures approach a new steady- 
state solution at time t = 0.3 s. (Fig. 11(d)), two linear 
conventional finite elements yield a fair approximation of the 
unknown nodal temperature but crudely approximate the temperature 
distribution. Both the nodeless parameter finite elements and the 
nodeless variable finite elements provide excellent prediction of 
the unknown nodal temperature and details of the nonuniform tempera- 


ture distribution. 



EXACT 

A- ▲ CONVENTIONAL FE. 



x/L 


(b) Comparative temperature distributions at t=0s, 


11. Conventional and nodeless variable finite element 
solutions for a rod with surface convection. 
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Even though the nodeless parameter approach employed in this 
problem yields excellent representation of the temperature distribu- 
tions at the beginning (t =* Os.) and near the end (t * 3.0 s.) 
of the response, the approach is unable to provide reasonable element 
temperature distributions during the response. As shown in Fig. 11(c) 
at time t * 0.01 s., temperatures obtained from the nodeless 
parameter finite elements are characterized by bumps within the 
elements. These unacceptable results are caused by using the steady- 
state element temperature distribution with the constant nodeless 
parameter for the transient analysis. Therefore, the nodeless 
parameter approach should not be employed for transient response 
predictions. Instead, the nodeless variable approach should be used 
since it gives accuracy superior to the linear conventional finite 
element throughout the response and predicts exact steady-state 
solutions, 

4.3.2 Transient Thermal Stresses in a Rod with Internal Heat 
Generation 

To further illustrate the use of the nodeless variable approach 
for one-dimensional transient problems and demonstrate additional 
benefits that can be achieved, an analysis of transient thermal 
stresses in a rod with internal heat generation is presented. 

A rod with constant cross-sectional area A and length L 

encased between fixed walls is shown in Fig. 12(a). Both ends of 

the rod have the specified temperatures at 311 K and 533 K at 

x = 0 and x * L, respectively. Initially, the rod is subjected 

3 

to a uniform internal heat generation rate Q = 353 kW/m and is 
in the thermal equilibrium. At time t = 0 + , 


the internal heat 
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Q 


1073 kW/m 3 



(a) ROD WITH INTERNAL HEAT GENERATION 


Fig. 12. Conventional and nodeless variable finite 
element solutions for a fixed end rod with 
internal heat generation. 
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3 

generation rate increases abruptly to 1073 kW/m and remains 
constant thereafter. The rod is modeled using: (1) 20 linear 

conventional finite elements, (2) two linear conventional finite 
elements, and (3) two nodeless variable finite elements. Comparative 
temperature distributions at time t « 0, 0.1 and 1.0 hr. are 

shown in Fig. 12(b). The figure shows that two nodeless variable 
finite elements have the same capability in predicting transient 
temperatures as 20 linear conventional finite elements. Two linear 
conventional finite elements underestimate the temperature distribu- 
tions with relatively large error throughout the transient response. 

In the structural analysis, three structural finite element 
nodels with the same discretizations as for the thermal finite 
element models are employed. Element temperatures obtained from 
the thermal finite element model are transferred directly to the 
structural finite element model for computation of displacements 
and stresses. For the quasi-static analysis, the structural response 
are computed at times corresponding to the transient thermal solu- 
tions obtained previously. At each time, the equivalent nodal 
thermal forces are computed using Eq. (3.57) and the element nodal 
displacements are computed from Eq. (3.55). Once the element nodal 
displacements are obtained, element displacement distributions and 
element thermal stresses are computed from Eqs. (3.53) and (3.58), 
respectively. 

Displacement distributions obtained from the three structural 
finite element models are shown in Fig. 12(c). The figure shows 
that two linear conventional finite elements are inadequate to 
represent the details of nonuniform displacement distributions. 
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(b) COMPARATIVE TEMPERATURE DISTRIBUTIONS 


Fig. 12. Conventional and nodeless variable finite 
element solutions for a fixed end rod with 
internal heat generation. 



R 

t = 0.1 HR 
t = 1.0 HR 


1 .25 
-2 
x 10* 

(c) COMPARATIVE DISPLACEMENT DISTRIBUTIONS 


Fig. 12. Conventional and nodeless variable finite 
element solutions for a fixed end rod with 
internal heat generation. 
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Similar to the thermal analysis, displacement distributions obtained 
from two nodeless variable finite elements and 20 linear conventional 
finite elements are in excellent agreement throughout the transient 
response. Comparative thermal stresses obtained from these finite 
element models at the times mentioned above are given in Table 5. 
Thermal stresses computed from two nodeless variable finite elements 
and 20 linear conventional finite elements are equal since temperature 
variations of these two finite element models coincided. Two linear 
conventional finite elements underestimate the thermal stresses and 
the error increases with time with a maximum of 10% at t = 1.0 hr. 

These two examples clearly demonstrate benefits of using the 
nodeless variable approach in one-dimensional transient thermal- 
structural problems. Further applications of the nodeless variable 
approach can be found in Ref. [28]. The use of the nodeless variable 
for improving temperature solutions in the transient thermal 
analysis directly improves accuracy of displacement and stress 
distributions in the structural analysis. The advantages of the 
nodeless variable approach for linear transient thermal-structural 
problems have been demonstrated in this chapter. The approach will 
be extended to nonlinear steady-state and transient thermal-structural 
analysis which includes radiation heat transfer in the next chapter. 
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Table 5 

Comparative Thermal Stresses for a Rod with 
Internal Heat Generation 


Time, t 
Hr. 


Stress, MPa 


2 Conventional 
Elements 

20 Conv. Elements 
2 Nodeless Variable 
Elements 

% Diff. 

0 

-507 

-531 

4.5% 

0.1 

-598 

-640 

6.5% 

1.0 

-652 

- 724 

10.0% 


i 


Chapter 5 


ONE-DIMENSIONAL THERMAL- STRUCTURAL FINITE ELEMENT ANALYSIS 
WITH RADIATION HEAT TRANSFER 

Due to their relatively low weight, high stiffness and ease of 
fabrication, trusses have high potential for use in space structures 
for solar collectors, antenna and space stations. Thermal analysis 
of these structures includes conduction heat transfer combined with 
significant radiation heat transfer. Radiation heat transfer intro- 
duces a strong nonlinearity in the energy equation being solved. 
Furthermore, a time dependent solution procedure is required for the 
analysis due to the changing orientation of the structure during the 
orbit. 

In this chapter, finite element solution procedures for one- 
dimensional transient thermal analysis with radiation heat transfer 
are presented. Three finite element types are formulated: (1) an 

isothermal element, (2) a linear conventional element, and (3) a 
nodeless variable element. Accuracy and efficiency of the finite 
elements are evaluated using two thermal-structural examples at the 
end of the chapter. 

5.1. Solution Procedures for One*-dimensional Transient Thermal 
Analysis with Radiation Heat Transfer 

In this section, transient thermal analysis for a one- 
dimensional finite element with radiation heat transfer is presented. 
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The radiation surface is assumed to be diffuse, gray and opaque 
which means the emitted radiation energy is uniformly distributed, 
independent of wave length and the material does not transmit 
radiation. For convenience all material thermal properties are 
assumed constant. 

For one-dimensional transient heat conduction in a rod with 
surface radiation, the governing differential equation for the 
temperature distribution T(x,t) can be derived using an energy 
balance on a small segment. With the assumptions mentioned above, 
the governing differential equation is 

2 

PcA |^ - kA + e crp s T 4 = a p q q’ r (5.1) 

0X 

where p is the density, c is the specific heat, A is the rod 

cross-sectional area, k is the material thermal conductivity, 

a is the Stefan-Boltzmann constant, e is the surface emissivity, 

a is the surface absorptivity, is the incident surface heating 

rate from distance directional sources per unit area, p and p 

s a 

are the cross-sectional perimeters for surface emitted energy and 
incident energy, respectively. 

Finite element equations corresponding to the governing differ- 
ential equation (5.1) can be derived using the method of weighted 
residuals as described in section 2.3. For this case, typical 
element equations have the form 

[C] {T} + [ [K c ] + [Kj] {T} * {Q c } + {Q r } (5.2) 

where [c] is the element capacitance matrix; [K^] and [K^_] 


are 
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the element conductance matrices corresponding to conduction and 
radiation, respectively, {Q^} is t * ie element vector of conduction 
heat flux across element boundary, and {Q r } is the element heat 
load vector due to incident radiation. These matrices are expressed 
in the form of integrals over the element length L as follows: 

EC3 - f 

0 

[K c ] - 5 

o 

L 

[K r ] {T} = J 
0 


L 

{Q r } = J ap q q r {N t } dx (5.4b) 

0 

where |_N^J denotes the element temperature interpolation functions. 

As shown in equation (5.3c), the conductance radiation matrix 
contains the element temperature within the integral. The element 
equations, Eq. (5.2), thus constitute a nonlinear set of equations. 
Since the time rate of change of the temperature vector {T} also 
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appears in the element equations, a transient nonlinear solution 
procedure is required for the analysis. 

Typical techniques for transient nonlinear solutions combine 
a linear transient solution method and a steady-state nonlinear 
solution method. The solution technique here uses a time-marching 
scheme where temperatures are computed at the middle of the time 
step, the Crank-Nicolson algorithm. At each time step, Newton-Raphson 
iteration is used to correct for nonlinearities. Further details of 
these methods including other solution algorithms can be found in 
the finite element text, Ref. [15]. 

Starting from the element equations, Eq. (5.2), the time- 
marching scheme is first applied by approximating the time rate of 
change of nodal temperatures as 


{T} 


At 


«T} 


n+1 


- (I> »> 


(5.5) 


where At is the time interval between the time step n and n+1 
such that t n+ ^ = t n + At; (T and {T> n+ ^ are the vectors of 
nodal temperatures at the time step n and n+1, respectively. 

Since the Crank-Nicolson algorithm computes temperature solutions at 
the middle of time steps, nodal temperatures at the middle of the 
step are approximated by 

(T } = j ( {T > n + (T} n+1 ) (5.6) 

where {T} denotes the vector of nodal temperatures at the middle 
of the step. From this equation, the vector of nodal temperatures 
at the step n+1 is 
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{T} = 2{T} - {T} (5.7) 

n-rl n 

By combining Eqs. (5-5) and (5.7), the time rate of change of the 
nodal temperature vector shown in Eq. (5.5) can be expressed in 
terms of nodal temperatures at the middle of the step and step n 
as 


{T} = -7 ({T} - {T} ) (5.8) 

At n 

Substituting Eq. (5.8) into Eq. (5.2), the element equations become, 

[77 [C] + [Kj + [Kl] {T} = (Q } + {Q } +77 [C] {T} (5.9) 

At ci. c r At n 

In Eq. (5.9), the vector of nodal temperatures (Tl^ that 
appears on the right-hand side is known from the previous step. 

Since the unknown nodal temperatures contained in the vector {T} 
are computed at the middle of time step, the heat load vectors must 
be evaluated at the same time. Once the unknown nodal temperature 
vector {T} is obtained, the nodal temperature vector { T } at 
the step n+1 can be computed from Eq. (5.7). 

The element equations obtained by applying the Crank-Nicolson 
algorithm shown in Eq. (5.9) are in the form of nonlinear algebraic 
equations 


[K(T)] (T> = {Q} 


(5.10) 


where 


[K(T)] {T} = [77 [C] + [Kj + [K r ] ] {T} 


(5.11a) 
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and 

(Q) - {Q c > + {Q r > +—■ [C] {T} n (5.11b) 

For any temperature vector (T) that is not an exact solution to 
equations shown in Eq. (5.10) above, an unbalanced nodal heat loads 
exist which can be written in the form of a vector as 

{*} - [K(T) ] {T} - (Q) (5.12a) 


or in tensor notations, 


r 
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j-1 


K. . T . - Q. 
ij 3 i 


(5.12b) 


To develop the Newton-Raphson method a Taylor series expansion with 
the first order-derivative accuracy is written as 


^({1}“) + c 




(5.13) 


A set of algebraic equations is obtained in the form 

[J]“ {AT} m+1 - {R} m (5.14) 

where the superscript m denotes the m t ‘ 1 iteration. The matrices 
[j] m and {R} m are the Jacobian matrix and the residual load vector, 
respectively, defined by 

3ip i 

J ij = “arT 


(5.15a) 
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R. 


1 



(5. 15b) 


At each iteration, the vector of nodal temperature increments 
{AT} m ~^ is computed using Eq. (5.14) and a new temperature vector 
is obtained from 


{T} m+1 = {T} m + {AT} m+1 (5.16) 

The iteration process is terminated when a convergence criteria 
(such as the maximum nodal temperature increment is less than a 
specified value) is met. For steady-state analysis, the equations 
shown in Eq. (5.2) do not contain the time rate of change of nodal 
temperatures, and only the Newton-Raphson iteration is required. 

5.2 Element Formulations 

In this section, three one-dimensional finite elements with 
surface radiation are formulated. Crank-Nicolson and Newton-Raphson 
methods described in the preceding section are employed for the 
transient and nonlinear solutions, respectively. 

5.2.1 Isothermal Element 

The isothermal element is a simple finite element suitable for 
problems with negligible conduction heat transfer. A uniform 
temperature variation is assumed along the element that varies only 
with time (Fig. 13(a)). This element is different from the finite 
elements mentioned in the previous chapters since element temperature 


is the only unknown for the isothermal element. Since the element 
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X 

(a) Isothermal element 




Fig. 13. One-dimensional element interpolation 
functions for nonlinear transient 
analysis with radiation. 
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neglects heat conduction, the governing differential equation shown 
in Eq. (5.1) becomes 

at 4 

pcA dt + £ ° P s T = a p q p r (5.17) 

where T denotes the element temperature which is a function only 
of time t. 

The Crank-Nicolson algorithm is applied to the above differential 
equation by first writing the rate of change of the element tempera- 
ture in the form of Eq. (5.8), 



(T - T n ) 


where T denotes the element temperature at the middle of the step. 
Substituting this equation into Eq. (5.17) yields a nonlinear 
algebraic equation in the form 

pcA + e ap s T 3 ) T = a p q q r + pcA T n (5.18) 

After the element temperature T at the middle of the step shown in 
the above equation is obtained, the element temperature at the end 
of the step is computed from 


T , - = 2T - T (5.19) 

n+1 n 

Next, the nonlinear algebraic equation shown in Eq. (5.18) is 
solved by applying the Newton-Raphson method. In this case, the 
unbalanced element heat load is given by (see Eq. (5.12)), 
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« ( -rz PcA + top T 3 ) T - a p q - ~ PcA T (5.20) 

v At *s r q n r At n 7 

Using Taylor series approximation, an algebraic equation is obtained 
in the form 


j“ AT™* 1 - R m (5.21) 

where J m and R m are the Jacobian and the residual load at the 
m th iteration defined by 

J® “ ff = Al pcA + 4 £ap s t3 (5.22a) 

R m * = -A pcA (T n - T) - ecr p g T 4 

(5.22b) 

+ a Pq p r 

At each iteration, the element temperature increment AT is computed 
from Eq. (5.21) and a new element temperature is obtained from 

T m+1 = T a + (5.23) 

After the convergence criterion is met, the element temperature 
shown in Eq. (5.23) is used in Eq. (5.19) to compute the element 
temperature at the end of the step. The use of the isothermal element 
does not require a set of simultaneous equations due to the assump- 
tion of negligible heat conduction as previously mentioned. The 
transient response of each element, therefore, can be computed 
separately. 
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The isothermal element is useful for modeling truss members 
where heat conduction is negligible in comparison with the incident 
heating and emitted radiation. Applications of the isothermal 
element for transient analysis of truss-type structures with surface 
radiation can be found in Refs. [29,30]. 


5.2.2 Conventional Element 

For the conventional element, a linear temperature variation 
is assumed between the two element nodes (Fig. 13(b)), 


T(x,t> = |l -f f J 1 


f = Ln t J m 


(5.24) 


[ T 2 (t) J 

where the unknown nodal temperatures T^(t) and ^(t) are a 
function of time t. With the conduction-radiation differential 
equation shown in Eq. (5.1), element equations can be derived as 
shown in Eq. (5.2). The vector of unbalance nodal heat loads shown 
in Eq. (5.12a) is written explicitly in the form 


( + > - [C](T} + [Kjm + [KJ(T} - (Q c ) - {Q r l 


- fil tCHT} 


n 


or 


{<M = -1 [c] { (T } - {T } } + [K ]{T> + [K Ht> 

At n c r 


- (Q c } - {Q r } 


(5.25) 
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For convenience nodal heat load vectors corresponding to each term 
on the right-hand side of the above equation are introduced to yield 


(<p) = {lp c } + {lp K ^} + “ {*Pq } 



(5.26) 


The Jacobian matrices and the residual heat load vector can now be 
formulated by using the definitions (see Eq. (5.15)) 



and 


R. = - 

i 1 


For example, the first term on the right-hand side of Eq. (5.26) is 
the heat load vector associated with the capacitance matrix. 


} = a7 [c]{{T> - {T} } 
At n 


L 

= At f PCA {N T } ^ n t J dx {{T} " {T} n } 
0 


With the linear element interpolation functions shown in Eq. (5.24), 
this term can be evaluated in closed form as 
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Using the definition of Jacobian 
corresponding Jacobian matrix is 


J. . 
i J 


i,j = 1.2 


the 
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(5,27) 


Similarly, the Jacobian matrix associated with the conductance conduc- 
tion matrix is obtained in the form 



(5.28) 


The third term on the right-hand side of Eq. (5.26) is the heat load 
vector associated with the radiation matrix. 


or 


L 

= [K r ] {T} = f EOp $ T 4 {N t } dx 

J o 


L 

ip . ■ 1 e ct p T 4 N . dx 
1 \ S X 

J o 


Therefore, the corresponding Jacobian is 


L 

9^i r 2 

J ij ' 3 t 7 * M EO p s 1 N l Hj <•* 
2 0 



0 


With the linear element interpolation functions shown in Eq. (5.24) 
this Jacobian matrix can be evaluated in closed form as 
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15 £0 Ps L 


10T^ + 6T^T 2 +3T ± t2 +l| 


2Tj +3 t£t 2 +31^ +2T| 


2T| +3T^T 2 +3X^2 +2T| 


+3T^T 2 +6X^2 +10X| 


(5.30) 


It can be seen that the Jacobian matrix associated with the 
radiation conductance matrix is strongly nonlinear since the unknown 
nodal temperatures contribute to all terms in the matrix. The matrix 
is sometimes [31] approximated by lumping these terms together 
similar to the lumped capacitance matrix given in Eq. (4.6). The 
lumped Jacobian matrix results in a much simpler form with zero off- 
diagonal terms, 

J Kr^ “ 2 eap s L 

« 

From Eqs . (5.15) and (5.26), the total residual load vector 
is 


(5.31) 


{R} - - + + * 


(5.32) 


For example, the residual load vector associated with the radiation 


conductance matrix is 
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{ R Kr } = = - [Kj. ] {T} 


L 

■ - f e a p s T 4 

0 


{N t } dx 


Using the linear element interpolation functions shown in Eq. 
this residual load vector can be evaluated in closed form. 


{R K> - - 30 £Gp s L i 


4 3 9 9 3 U 

5T 1 + 4X^2 + 3T 1 T2 + 21^^ + 

+ 21^X2 + 3 T*t| + + 5 T^ 


(5.24), 


V (5.33) 


J 


After all Jacobian matrices and residual heat load vectors are 
computed, a final set of algebraic equations is obtained in the form 


[J] m { AT} m+1 = {R} m 


(5.34a) 


where the superscript m denotes the mth iteration and, 


[J] - [J c ] + [J K ] + [J K ] (5.34b) 

W - (R Cc > + ( r K c > + ( R K r > + ( r Q c > + {®Q r > (5.34c) 


The solution of the temperature vector at successive times proceeds 
as previously discussed for the isothermal element. 

As shown in Eq. (5.34a), the transient and nonlinear solution 
procedures lead to a set of algebraic equations. The Jacobian 
matrices and the residual heat load vectors shown in Eqs. (5.34b-c) 
are thus necessary for the analysis. With the linear element inter- 
polation functions shown in Eq. (5.24), these matrices can be 



115 


evaluated in closed form and are given as computer subroutines in 
Appendix D. 


5.2.3 Nodeless Variable Element 

In the preceding chapter, the nodeless variable approach was 
introduced for improvement of temperature solutions in one-dimensional 
linear transient analysis. The basic idea of the approach is the use 
of the steady-state element temperature interpolation function derived 
from the solution of a given ordinary differential equation. An 
element nodeless variable is employed so that exact steady-state 
solutions are obtained at the beginning and at the end of the 
transient with realistic temperature distributions prediced throughout 
the response. For one-dimensional conduction-radiation heat transfer, 
it is not possible to obtain a closed form solution to the governing 
differential equation, Eq. (5,1). However, the nodeless variable 
approach is still useful for the analysis to provide improved 
temperature solutions for the thermal element while maintaining the 
same discretization as the two node structural element. The element 
temperature distribution with a nodeless variable is written in the 
form, 


T(x,t) = lN 0 (x) N^x) N 2 (x)J 




T 0 (t) 

T l (t) [ 
T 2 (t) 


Ln t J m 


(5.35) 


where N^(x) is the nodeless variable interpolation function, 

N_^(x), i =1,2 are typical element interpolation functions; Tg(t) 
is the nodeless variable, and T_^(t), i =1,2 are the nodal tempera- 


tures . 
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As mentioned earlier, the nodeless variable interpolation func- 
tion Nq must vanish at nodes in order to preserve continuity of 
temperature between elements. There are wide choices for selecting 
the nodeless variable interpolation function to meet this requirement. 
The simplest function is in the form of polynomials with one order 
higher than the linear element interpolation functions used in the 
conventional finite element, 

N 0 (3 ° "§(!“?) (5.36a) 

N 1 (x) = 1 “ § (5.36b) 

N 2 (x) = f (5.36c) 

With these element interpolation functions, the element temperature 
distribution, Eq. (5.35), results in a parabolic distribution over 
the element length as illustrated in Fig. 13(c). 

The use of the nodeless variable element for transient heat 
conduction with surface radiation follows the same procedure 
described for the linear conventional element. Typical element 
equations derived from the method of weighted residuals shown in 
Eq. (5.2) contain three unknowns. These element unknowns are the 
nodeless variable T^ and two nodal temperatures T^ and T^. 

For transient solutions, the Crank-Nicolson algorithm is applied, 
and a set of nonlinear algebraic equations is obtained. Next the 
Newton-Raphson method is used and a new form of simultaneous 
Eq. (5.34a) is obtained where the Jacobian matrices and the residual 
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load vectors are defined by Eqs. (3.54b-c), respectively. For example, 
the Jacobian matrix contributed from the radiation conductance matrix 
has the form 


L 

[J K ] - 4 f e a p g T 3 {N x > |nJ dx 
r J„ 


Using the element nodeless variable temperature distribution and 
their element interpolation functions shown in Eqs. (5.35) and (5.36), 
respectively, this Jacobian matrix is written explicitly as, 




- 4 ( P s (¥o + S 1 T 1 + B 2 T 2 ): 

J r\ 


l L N 0 N i n 2 J dx 


Due to the complexity of the Jacobian matrix as shown above and 
other matrices that appear in Eq. (5.34), the computer-based symbolic 
manipulation language MACSYMA [32] was used to evaluate the matrices 
in closed form. Results of the Jacobian matrices and the residual 
load vectors are provided in the form of computer subroutines in 
Appendix D. 

After the Jacobian matrices and the residual load vectors are 
computed, typical element equations shown in Eq. (5.34a) can be 
written in the form, 
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These element equations contain unknowns in the increments of the 
nodelss variable and the nodal temperatures, i.e. one more unknown 
than those obtained from the linear conventional element* Once these 
unknowns are obtained, new values of the nodal temperatures and the 
nodeless variable are computed using Eq. (5*16). After the iteration 
process is terminated, the nodal temperatures and the nodeless 
variable at the end of time step are computed from Eq. (5.7). 

Finally, the temperature distribution within the element is computed 
by using the element nodeless variable interpolation functions shown 
in Eq. (5.35). 

It should be noted that the nodeless variable interpolation 
functions, Eq . (5.36), introduced in this section are applicable 
when other heat transfer modes (such as surface convection) are 
included in the analysis. The element temperature distribution in 
the parabolic form can provide a more realistic temperature distribu- 
tion than the linear conventional element. This type of the nodeless 
variable interpolation functions suggests that the nodeless variable 
approach can be generalized to other finite element types. To 
investigate this possibility, a two-dimensional nodeless variable 
thermal element is developed in the next chapter. 

5.3 Applications 

The effectiveness of the nodeless variable finite element 
described in this chapter is demonstrated for two examples of conduc- 
tion and radiation heat transfer. The linear conventional finite 
element described in section 5.2.2 is used in these two examples for 
comparison of solution accuracy. Temperatures computed from the 


119 


nodeless variable and linear conventional finite elements are used 
in the structural analysis for computation of displacements and 
thermal stresses. 

5.3.1 T hermal Stress in a Rod with Surface Radiation 

A rod with constant cross-section area A and length L 
encased between fixed walls is shown in Fig. 14(a). The rod has 
specified end temperatures at 311 K and 533 K at x « 0 and 
x - L, respectively, and is cooled along the surface by radiation 
to zero medium temperature. The rod is modeled using (1) 20 conven- 
tional elements with consistent Jacobian matrices (see Eq. 5.30)), 

(2) two conventional elements with consistent Jacobian matrices, 

(3) two conventional elements with lumped Jacobian matrices (see 
Eq. (5.31)), and (4) two nodeless variable elements. The terms 
consistent and lumped refer to the formulation of the Jacobian matrix 
contributed by the radiation conductance matrix described in section 
5 . 2 . 2 . 

Temperature distributions computed from these four finite 
element models are compared as shown in Fig. 14(b). The figure shows 
that two nodeless variable elements have the same capability in 
predicting the unknown nodal temperature (at x/L = 0.5) and element 
temperature distributions as 20 conventional finite elements. Two 
conventional finite elements with consistent formulation underesti- 
mate the unknown nodal temperature and crudely approximate temperature 
distribution. Two conventional finite elements with lumped formula- 
tion overestimate both the unknown nodal temperature and element 
temperature distributions with relatively large error. 


RADIATING HEAT 



(a) ROD WITH SURFACE RADIATION 


Fig, 14. Conventional and nodeless variable finite 
element solutions for a fixed end rod 
radiating to space. 
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X/L 

<b) COMPARATIVE TEMPERATURE DISTRIBUTIONS 


Fig. 14. Conventional nodeless variable finite element 
solutions for a fixed rod radiating to space. 
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In the structural analysis, four structural elements with the 
same discretizations as for the thermal models are employed. Element 
temperatures obtained from the thermal model are transferred directly 
to the structural finite element model for computation of displace- 
ments and stresses. The conventional structural finite elements 
employ linear element displacement distributions as used in typical 
finite element programs. The structural finite element for the 
nodeless variable thermal element uses the exact displacement 
distribution, Eq. (3.53), derived based upon the parabolic element 
temperature distribution shown in Eq. (5.35). Displacement distribu- 
tions obtained from these structural finite element models are 
compared as shown in Fig. 14(c). The figure shows that two conven- 
tional finite elements are inadequate to represent the nonuniform of 
displacement distribution. In addition, two conventional finite 
elements with consistent and lumped formulations overestimate the 
thermal stress (not shown) by 12 and 23 percent, respectively. 
Displacement distributions obtained from two nodeless variable 
finite elements and 20 conventional finite elements are in excellent 
agreement where the difference in the thermal stresses is negligible 
(less than 0.05 percent). 

5.3.2 Thermal Analysis and Structural Response of a Space 
Truss Module 

A three member orbiting truss module shown in Fig. 15(a) is 
used to demonstrate the efficiency of the nodeless variable finite 
element. A typical truss member receives incident heating which is 
a combination of: (1) solar heating, (2) earth emitting heating, 

and (3) earth reflected solar heating. With the open-truss type 





(a) Orbiting truss space structure. 


Fig. 15. Thermal analysis and structural response of a space truss module. 
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structure as shown in the figure, member to member radiation exchanges 
are relatively small and are neglected. A geosynchronous orbit 
(period of 24 hr.) is employed where solar heating is large compared 
to the earth emitted heating. During the orbit, incident heating 
normal to a typical truss member varies continuously due to the 
changing orientation of the member. As the orbiting truss enters 
and leaves the earth’s shadow, the incident heating changes rapidly. 
Member temperatures and structural deformations thus depend strongly 
on the time-dependent incident heating. 

To demonstrate the use of the conventional and the nodeless 
variable finite elements formulated in the preceding section, the 
truss module with properties of aluminum is considered. Four finite 
element models are employed where each truss member is represented 
by: (1) 10 conventional elements with consistent formulation, (2) one 

conventional element with consistent formulation, (3) one conventional 
element with lumped formulation, and (4) one nodeless variable 
element. Temperature distributions computed from these four finite 
element models at a typical orbital position are shown in Fig. 15(b). 
The figure shows that the nodeless variable finite element model 
provides excellent prediction of the nodal temperatures and very 
good element temperature distributions compared to the refined 
conventional finite element model. The conventional finite elements 
with consistent formulation tend to average the nonuniform tempera- 
ture distributions and thus cannot provide accurate nodal temperatures. 
The conventional finite elements with lumped formulation predict 
nodal temperatures very well but yield large errors for member 
interior temperatures. Comparative temperature distributions of 
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(b) Comparative temperature distribution of a three member 
orbiting truss at 6 = 60 degrees. 


Fig. 15. Thermal analysis and structural response of a space truss module. 
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finite element models at other orbital positions also show a similar 
trend; the nodeless variable finite elements predict nodal tempera- 
tures and member temperature distributions very accurately compared 
to the refined conventional finite elements. 

Temperature obtained from the four thermal finite element models 
for a complete orbit are transferred to the structural finite element 
models for computation of displacement histories. The quasi-static 
analysis described in section 2.4 is employed for the computation 
of the unknown nodal displacements. Fig. 15(c) shows a comparison 
of typical member elongation histories computed from the finite 
element models during the orbit. Since the temperature distributions 
obtained from the nodeless variable finite element model and the 
refined conventional finite element model are in very good agreement, 
member elongation histories predicted by these two finite element 
models almost coincide (maximum difference of 1 percent) . Conven- 
tional finite element models with consistent and lumped formulations 
yield errors for member elongation up to 29 and 44 percent, 
respectively. Such large errors result from the incapability of the 
conventional finite element to provide realistic member temperature 
distributions. Since the conventional finite element with consistent 
formulation trends to average the member temperature as previously 
mentioned, accuracy of the member elongation history obtained is thus 
higher than the conventional finite element with lumped formulation. 

These two examples clearly demonstrate the benefits of using 
the nodeless variable finite elements in one-dimensional radiation- 
conduction problems that are characterized by nonuniform temperature 
distributions. The elements predict member temperatures accurately 
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(c) Comparative displacements of a three member orbiting truss. 


Fig. 15. Thermal analysis and structural response of a space truss model. 


128 


129 


and are compatible with two-node structural elements to permit an 
integrated thermal-structural analysis. Additional applications of 
the nodeless variable finite element for one-dimensional thermal 
problems with conduction and radiation heat transfer appear in 
Ref. [33]. 



Chapter 6 


TWO-DIMENSIONAL NODELESS VARIABLE FINITE ELEMENTS 

In the two preceding chapters the nodeless variable approach 
was applied to one-dimensional linear thermal-structural analysis 
and to nonlinear radiation heat transfer. The unique feature of the 
approach is the use of an additional nodeless variable for a thermal 
finite element. Improvement of solution accuracy is achieved while 
the same discretization is employed for both thermal and structural 
finite element models. 

In this chapter the approach is extended for development of 
two-dimensional nodeless variable finite elements. Restrictions for 
developing these finite elements are first discussed. Two nodeless 
variable finite elements and their interpolation functions are 
presented. Then the use of the nodeless finite elements for linear 
thermal-structural analysis is described. Efficiency of the nodeless 
variable finite elements is evaluated by comparison with the 
conventional bilinear four-node finite element and exact solutions 
in examples at the end of the chapter. 

For simplicity in understanding characteristics of the two- 
dimensional nodeless variable finite elements, a brief description 
of a conventional bilinear four-node thermal finite element is 
first given. The element temperature distribution for a bilinear 
four-node element is expressed in the form. 
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t - Ln x n 2 n 3 n 4 J 


= Ln t J {t} 


( 6 . 1 ) 


where N^, i =1,4 are the element interpolation functions which 
are a function of spatial coordinates in two-dimensions, and T\, 
i =1,4 are the time dependent nodal temperatures. 

Fig. 16(a) shows a conventional four-node element with a general 
quadrilateral shape. As described in section 2.3, typical finite 
element matrices are in the form of integrals over the element 
volume or along the element boundary. Such element matrices for 
a quadrilateral shape are difficult to evaluate. To simplify the 
element integrations, the quadrilateral element in the Cartesian 
coordinate system (x,y) is transformed to a natural coordinate 
system (£,n) as shown in Fig. 16(b). The two coordinate systems 
are related by 


4 

x = Z N . (£,n) x. = [Nj {x} (6.2a) 

i=l 

4 

y = Z N . (C,n) y. = L N J (y> (6.2b) 

i=l 1 


where N_^, i =1,4 are the element shape functions defined by. 


»1 - jd-O (l-n) 
n 3 = i(l+g) (1+n) 


n 2 = i(i+o (l-n) 
n 4 = ^(1-5) (l+n) 


(6.3a) 
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or in compact form, 

N i " \ (1 + 5C i )(1 + nn i ) i -1,4 (6.3b) 

where and i =1,4 are the nodal coordinates in the natural 

coordinate system. For example, 5^ = = -1, ^ = n 2 = 

etc. 

When the shape functions shown in Eq. (6.3) are used as 

the element temperature interpolation functions in Eq. (6.1), this 
conventional element is called an isoparametric quadrilateral element 
because the same interpolation functions are used to interpolate 
temperature and spatial coordinates. 

Note that an element temperature interpolation function shown 
in Eq. (6.3) has a value of unity at the node to which it pertains 
and a value of zero at the other nodes. Along the element edge 
(£ = ±1, n = ±1), these element interpolation functions are linear. 
Therefore, the temperature distribution along a typical element 
edge varies linearly where the magnitude depends on the temperatures 
of the two corner nodes located at that edge. When elements are 
connected, the conventional quadrilateral element preserves 
continuity of temperature along the element interfaces. The conti- 
nuity of the element interface temperatures is a basic requirement 
to assure convergence of the temperature solution as element size 
decreases. This continuity requirement must be met when a new 
thermal finite element is constructed. Further details of require- 
ments for a typical finite element to meet convergence criteria can 
be found in Ref. [15] . 
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6.1 Two-Dimensional Nodeless Variable Thermal 
Finite Elements 

In several thermal-structural applications, a more detailed 
finite element thermal model is required than the finite element 
structural model. To maintain the same discretization for thermal 
and structural models, new thermal finite elements are required. In 
this section, two type of two-dimensional nodeless variable thermal 
finite elements for improved temperature solutions are presented. 
These elements predict more realistic temperature distributions than 
the conventional finite element previously described. The basic 
objectives for developing the new finite elements are: (1) the 

elements should provide a nonlinear temperature distribution but 
maintain four element nodes to be congruent with the four node 
structural element, and (2) compatibility of temperature along 
element interfaces must be preserved. The nodeless variable concept 
previously described for one-dimensional element is extended to 
two-dimensions to meet these objectives. 

6.1.1 "Bubble" Nodeless Variable 

One approach [34] for constructing nodeless variable finite 
elements is to add a "bubble" function which vanishes along the 
element boundaries. The element temperature distribution is written 
in the form, 

1 
2 

3 

4 
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where is the nodeless interpolation (bubble) function defined 

by 


n q = (1 - c 2 )d - n 2 ) 


(6.5) 


and Tq is the nodeless variable. 

Along the element boundary (5 = ±1, n = ±1), the bubble 
function, Eq. (6.5), is identically zero. Therefore, the element 
boundary temperature reduces to a linear variation as for the 
conventional finite element and continuity of temperature along 
element interfaces is preserved. Within the bubble nodeless variable 
element, the temperature distribution is a combination of the 
conventional element temperature distribution and a bubble function 
where its magnitude is measured by the nodeless varaible Tq. The 
combination thus permits a quadratic temperature distribution over 
the element. 

It should be noted that even though the bubble nodeless 
variable finite element can provide a quadratic temperature distribu- 
tion within the element, the temperature along the element boundary 
is linear. To achieve further improvement of the temperature solu- 
tion, the temperature distribution should vary nonlinearly along the 
element boundary. With the idea of the bubble function, a nodeless 
variable finite element with this behavior can be constructed. This 
type of nodeless variable finite element is presented in the next 


section. 
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6.1.2 Boundary Nodeless Variable 

In order to establish a nonlinear temperature distribution along 
the four element edges as well as within the element interior, the 
following four nodeless variable interpolation functions (see Fig. 17) 
are employed 


N 5 = \ (1 “ ?2)(1 ' n) (6.6a) 

N 6 = y (1 + 5) (1 - n 2 ) (6.6b) 

N? - \ (1 - C 2 )(l + n) (6.6c) 

Ng = (1 - 5)d - n 2 ) (6.6d) 


where each interpolation function varies quadratically along one 

edge and vanishes on the other edges. For example, the nodeless 

2 

variable interpolation function varies as 1-5 along the 

edge n = -1 and is identically zero on the other three edges. 

As mentioned earlier, continuity of the temperature along the 
element interfaces must be assured for convergence of the solution. 
This restriction can be met by providing a nodeless variable for 
each element edge. With a nodeless variable for each element edge, 
element interpolation functions for a quadrilateral element can be 
written in the form, 
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(6.7) 



n 7 = |(l-r)(l+n) 


Ng = |(l-C)(l-n ) 


Fig. 17. Nodeless variable interpolation functions for 
two-dimensional quadrilateral finite element. 
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where T\, i =1,4 and i =5,8 are the nodal temperatures and the 
nodeless variables, respectively. Element interpolation functions 
N., i =1,4 are the same as for the conventional bilinear four node 

i 5 

element given in Eq. (6.3), and N^, i=5,8 are the nodeless variable 
interpolation functions given in Eq. (6.6). 

The combination of the conventional and nodeless variable inter- 
polation functions, Eq. (6.7), provides a quadratic temperature 
distribution over the element but with only four element nodes. 
Interelement compatibility is preserved since adjacent elements have 
a common nodeless variable on adjoining edges. The magnitude of the 
nonlinear variation on an element edge is measured by the correspond- 
ing nodeless variable. Temperature distributions for the conventional 
bilinear element and the nodeless variable element are compared in 
Fig. 18. 


6.2 Nodeless Variable Finite Element Formulation 
for Thermal-Structural Analysis 

In this section, the thermal finite element formulation for two- 
dimensional linear transient analysis is described. The formulation 
is valid for both the conventional element and the nodeless variable 
element. A four node structural element which will be used in 
junction with the thermal element for computation of thermal stresses 
is also presented. 

6.2.1 Linear Thermal Analysis 

In two-dimensional transient heat conduction, the governing 
differential equation for the temperature distribution T(x,y,t) 
may be expressed in the form of 
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Fig. 18. Two-dimensional element interpolation 
functions. 
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( 6 . 8 ) 


where k and k are the thermal conductivities in x and y 
x y 

directions, respectively, Q is the internal heat generation rate 
per unit volume, p is the density, and c is the specific heat. 

To derive the element equations and element matrices, the method 
of weighted residuals (see section 2.3) is applied to the governing 
differential equation (6.8). With the boundary conditions of 
specified temperatures, surface heating and surface convection as 
shown in Eqs. (2.5a-c), typical element equations have the form 

[C] {f} + [K c + K^] (T> = { Q c > + {Qq} + {Q q } + {Q h > (6.9) 

where [c] is the element capacitance matrix; [K c ] and [k^3 are 
element conductance matrices corresponding to conduction and convec- 
tion, respectively. These matrices are expressed in the form of 
integrals over the surface area A of an element with the thickness 
t as follows: 


[c] = t | pc {N t } L n t J dx dy 


(6.10a) 


[Kj = t J [B t ] T [k] [B T ] dx dy 


(6.10b) 


[K^] = _[ h{N x } [N X J dx dy 


(6.10c) 


where [B^] denotes the temperature gradient interpolation matrix, 
and h is the convection coefficient. The right-hand side of the 



discretized equation (6.9) contains heat load vectors due to 
specified nodal temperatures, internal heat generation, surface 
heating, and surface convection. These vectors are defined by 
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{Q c > 

V 





(q • h) {N t } ds 

1 

Q{N t } dx dy 
q(N^} dx dy 


J h T a (N t ) dx dy 


(6.11a) 


(6.11b) 


(6.11c) 


(6. lid) 


where q is the vector of conduction heat flux across boundary 
that is required to maintain the specified nodal temperatures, q is 
the surface heating rate per unit area, and is the convective 

medium temperature. 

As mentioned earlier, a typical quadrilateral element in 
Cartesian coordinates (x,y) is transformed to the natural coordi- 
nates (£,n) to perform the element matrix integration. In 
computation of the conduction conductance matrix (Eq. (6.10b)), for 
example, the chain rule Is first applied to relate the temperature 
gradients in both coordinate systems. 
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( 6 . 12 ) 
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Using the coordinate trasnformation shown in Eq. (6.2 
relations become 
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where [j] is the Jacobian matrix defined by 


[J] = 


4 3N. 


4 3N 
E 


i=l 


3n X i 


4 3N. 

£ ^ y ‘ 


4 3N. 

i-1 9n 


Substituting the element temperature, Eq. (6.1) or (6 
right-hand side of Eq. (6.13) yields 


3T 

3x 

3T 

3 y 


= [b t (5,h)] 


T 

r 

v- J 


where r is the number of the element unknowns; r = 


, the above 


(6.13) 


(6.14) 


.7), into the 


(6.15) 


4 and 8 for 


the conventional bilinear element and the nodeless variable 
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element, respectively. The temperature gradient interpolation matrix 
in the above equation is given by 


3N X 3N 2 

~W IT 


u t (e,ti)] - ur 1 


8^ 3Nj 

3n 3n 



3N 
r 

an 


(6.16) 


Using dx dy = |j| d£ dn where |j| is the determinant of [j] , 
the conduction conductance matrix terms of the natural coordinates 
is 


[K c ] = t 


[B T (?,n)j 


[k] [B T (5,n>] |j| dn (6.17) 


-1 -1 


Next, the coefficients in the conduction conductance matrix 
are computed by numerical integration; the Lagendre-Gauss method 
is used where the above conduction conductance matrix is written in 
the form, 

NG NG 

[K c ] - ^ ^ w ± w. [B T (5 i ,n j )] i [k] [B T (5 i ,n J )] ( 6 . 18 ) 

where W^, W_. denote Gauss weight factors, rij denote gauss 

integration points and NG is the number of Gauss points in each 
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coordinate direction. Gauss weight factors a,nd Gauss integration 
point coordinates can be found in Ref. [ 15] - 

Other element matrices shown in Eqs . (6.10 -6.11) can be 
formulated in the same manner. For example, the conductance 
matrix and the heat load vector associated with surface convection 
are expressed as 


NG NG 

[k. ] = h r z 

i=l j=l 


w. w. 

1 3 


[N T (c i ,n J .) 3 


[N T a., 


n. )] 




(6.19a) 


NG NG 
{Q } = Z Z 

i=l j=l 


w. w. 

1 3 


[N T (5 i9 rij)] T |j(5 i ,n^)| 


(6.19b) 


In performing the numerical integration, the accuracy of the 
matrices depends on the number of Gauss points used. In general, 
the use of n Gauss points provides exact integration when the 
integrand contains polynomials of order up to 2n -1. For the 
conventional bilinear finite element, two Gauss points (NG = 2) 
in each coordinate direction are normally used. Since the nodeless 
temperature interpolation functions contain higher order of 
polynomials than those for the conventional bilinear element, 
more Gauss points should be used. For the linear thermal analysis 
presented herein, three Gauss points in each coordinate direction 


145 


was found by numerical tests to be appropriate for accurate node- 
less variable element matrices. 

After element matrices are computed, typical element equations 
can be written in the form 


[c] {T} + [K] {T} = {Q} (6.20) 

The conventional bilinear element has four nodal temperatures as 
the unknowns, thus the above element equations contain four 
unknowns. The nodeless variable element has four nodal tempera- 
tures and four nodeless variables as the element unknowns, therefore, 
the element equations contain eight unknowns. In transient 
analysis, these eight equations must be solved simultaneously 
similar to the one-dimensional nodeless variable described in the 
preceding chapter. In steady-state analysis, the four nodeless 
variable unknowns can be eliminated from the element equations 
using the matrix condensation technique [35]. The final number 
of element equations thus reduces to be the same as of the conven- 
tional bilinear element. 

6.2.2 Structural Element 

In this section, the congruent structural element is briefly 
described. The element contains four nodes and permits the same 
discretization with the conventional and nodeless variable thermal 
elements described in the preceding sections. The element stiff- 
ness matrix is the same as used in conventional four node structural 
elements. However, the improved element temperature distributions. 
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Eq. (6.7), are incorporately consistently in the thermal force 
vector computation to yield an integrated thermal structural 
element. 

The structural element at each node has two in-plane displace- 
ment unknowns u and v which may vary with the element local 
coordinates x, y and time t. Element displacement distributions 
are assumed in the form (see Eq. (2.23)), 
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o n 3 
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O 

CM 

2 ; 

n 3 0 n 4 


u„ 


u„ 


= [N g ] {&} 


( 6 . 21 ) 


where N^, i =1,4 are the element displacement interpolation 
functions which have the same form as for the conventional finite 
element temperature interpolation functions shown in Eq. (6.3). 

For the quasi-static analysis, typical element equations shown 
in Eq. (2.26) reduce to 


[K g ] { 6 } - (F t > 


where [K 1 
s 


is the element stiffness matrix, and {F^} is the 


( 6 . 22 ) 



equivalent nodal thermal load vector. These matrices are expressed 
in the form of integrals over the element volume V as 
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[D] [B ] dv 
s 


{F x > = | [B g ] T [D] {a} (T (x,y , t) - T. ef ) dV 

where [B g ] is the strain-displacement interpolation matrix 
from the strain-displacement relations. 


3u 

3x 


xy 


3v 

3y 

3u 3 v 
3y 3x 


[B s ] {&} 


[d] is the elasticity matrix defined by (plane stress) , 

[D] = — ^2 
1-v 


where v is Poisson's ratio. The vector {a} contains the 
expansion coefficients given by (plane stress) 



{a} = < a 


(6.23a) 

(6.23b) 

obtained 

(6.24) 

(6.25) 
thermal 

(6.26) 


0 
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T(x,y,t) is the element temperature computed using the conventional 

or the nodeless variable thermal element and T is the reference 

ret 

temperature for zero stress. The elasticity matrix [D] and the 

vector of thermal expansion coefficients shown in the above equations 

2 

can be used for plane strain by replacing E/(l-v ) for E, v(l-v) 
for v, and (l+v)a for ct. 

Similar to the quadrilateral thermal element, numerical integra- 
tion is required for computing the element matrices. Using the 
Lagendre-Gauss method, the element stiffness matrix and the equivalent 
thermal load vector shown in Eqs. (6.23a-b) are written in the form, 


NG NG 

[k ] = t z 


z w w [B (£ n >r [d] [B (5 .,n.)] |J(5.,n ) 
i=l j=l J J J J 


(6.27a) 


NG NG 


{F x } = t 


W i W j [B s U i’ n j )] [ °] {a} (T( Si^j) -T ref ) |j(C.,n.) 


(6.27) 

where T(5^,rij) is the temperature at the element Gauss integration 
point 5^. and rij • 

Unlike the thermal finite element previously described, the nodal 
displacement unknowns of the structural element are the vector 
quantities. Transformation of the element matrices from the local 
coordinates (x,y) to the global coordinates (X,Y,Z) is required. 

In three-dimensions, the element stiffness matrix becomes a 12 by 12 
matrix and similarly with the nodal force vector. Thus the element 
equations contain a total of 12 equations with 12 nodal displacement 
unkonwns in the global coordinates. After the global element matrices 
are assembled and the nodal displacements are computed, element nodal 
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displacements in the local coordinates can be obtained. Then the 
element stresses can be computed from 


r 


< 



[D] 


[B ] {6} - { ct } 

S 


(T(x,y,t) 



(6.28) 


6.3 Applications 

To illustrate the performance of the two-dimensional nodeless 
finite element presented in section 6.1.2, two examples are analyzed: 

(1) a rectangular plate with surface convection, and (2) a simplified 
wing section with aerodynamic heating. In each example, benefits 

of the nodeless variable finite element are demonstrated by comparison 
with results from conventional finite element and analytical solutions. 

6.3.1 Steady-State Heat Conduction in a Plate with Surface 
Convection 

A rectangular plate (Fig. 19(a)) has a specified temperature 
Tq along the boundaries. The plate is cooled by surface convection 
to a zero medium temperature, T^ « 0. Using symmetry, a quarter of 
the plate is first modeled by: (1) one conventional element, and 

(2) one nodeless variable element. 

Fig. 19(b) shows -the comparative temperature distributions at 
y = b/2 for an analytical solution [27] , the conventional element 
and the nodeless variable element solutions. For these models, the 
conventional element gives a relatively high error compared to the 
nodeless variable element. The largest error for both finite element 
models occurs at the center of the plate (16% and 3% for the 
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ANALYTICAL 

▲ ▲ CONVENTIONAL F.E., I ELEMENT 

O- — -O NOOELESS VARIABLE F.E., I ELEMENT 



(b) COMPARATIVE TEMPERATURE DISTRIBUTIONS ALONG y * 


Fig. 19. Conventional and nodeless variable finite element 
solutions for a plate with surface convection. 
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conventional element and the nodeless variable element, respectively). 
At the center of the plate, both elements show a discontinuity of 
conduction heat flux indicating a need for mesh refinement. Next, 
the plate is modeled by using four finite elements shown by the dotted 
lines in Fig. 19(a); comparative temperature distributions are shown 
in Fig. 19(c). Four conventional elements provide a fair estimate 
of the temperature variation, but four nodeless variable elements 
yield excellent predictions for both nodal and element temperatures. 
Comparisons of temperatures at other sections of the plate (not shown) 
demonstrate that four nodeless variable elements provide excellent 
agreement with the analytical solution for the entire plate. 

6.3.2 Simplified Wing Section with Aerodynamic Heating 

To demonstrate the usefulness of the two-dimensional nodeless 
variable elements in aerospace thermal-structural analysis, a simpli- 
fied wing section is analyzed (Fig. 20(a)). Top and bottom skins 
of the wing section are connected by three corrugated spars and are 
subjected to symmetrical, nonuniform time-dependent aerodynamic 
heating. 

Three finite element models are employed to computed temperatures. 
For a unit depth in the spanwise direction, the first model consists 
of seven conventional elements; two elements each for the top and 
bottom skins and one element for each spar. The second model is 
identical to the first model except nodeless variable elements are 
used. The third model uses a refined mesh (not shown) with 35 conven- 


tional elements. 
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(c) COMPARATIVE TEMPERATURE DISTRIBUTIONS ALONG y * -|- . 


Fig. 19. Conventional and nodeless variable finite element 
solutions for a plate with surface convection. 
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(a) SIMPLIFIED WING SECTION WITH AERODYNAMIC HEATING 


Fig. 19. Conventional and nodeless variable finite element 
solutions for a simplified wing section with 
aerodynamic heating. 



r 
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Comparative skin temperature distributions at t = 150 s. are 
shown in Fig- 20(b); the number of elements cited is for the skin 
only. The nodeless variable finite element model predicts a realistic 
temperature distribution and gives very good agreement with the 
result from the refined conventional finite element model. The crude 
conventional finite element model underestimates the average skin 
temperature and is unable to provide details of the nonuniform 
temperature distribution. 

In computation of the skin thermal stress, classical beam 
theory [16] is employed for comparison with two finite element stress 
analyses- Detailed temperature distributions from the refined 
conventional finite element thermal model are used to compute the 
stress a x from beam theory. Temperature distributions from the 
crude conventional thermal finite element model and the nodeless 
variable thermal finite element model are transferred to a structural 
finite element model with the same discretization for the stress 
computations. Comparative stress distributions at t = 150 s. are 
presented in Fig. 20(c). The advantage of using the improved 
temperature distributions from the nodeless variable finite element 
model in computing stresses is clearly demonstrated. These stress 
distributions are in excellent agreement with the result from beam 
theory with both the critical stress and its location accurately 
predicted. Using the temperature distribution from the crude 
conventional finite element model yields significant errors in the 
stress distribution and is unaccetable for this problem. 

These two examples clearly demonstrate the benefits of the 


two-dimensional nodeless variable finite element that can be obtained 
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(b) COMPARATIVE TEMPERATURE DISTRIBUTIONS 
ALONG CHORDWISE DIRECTION, t = 150 s. 


Fig. 20. Conventional and nodeless variable finite element 
solutions for a simplified wing section with 
aerodynamic heating. 



(y,l50), MPa 
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(c) COMPARATIVE THERMAL STRESS DISTRIBUTIONS 
ALONG CHORDWISE OIRECTION, t = 150 s. 


Fig. 20. Conventional and nodeless variable finite element 
solutions for a simplified wing section with 
aerodynamic heating. 
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for thermal-structural analysis. Additional applications, a summary 
of the nodeless variable approach and the thermal-structural finite 
element formulation presented in this chapter appear in Ref. [36]. 



Chapter 7 


CONCLUDING REMARKS 

An integrated approach for improved thermal-structural finite 
element analysis is presented. The approach was motivated by aero- 
space applications to improve thermal-structural finite element 
analysis capabilities. An important goal is to eliminate the 
incompatibility between thermal-structural analyses where a more 
detailed finite element model is required for the thermal analysis 
than for the structural analysis. The integrated approach is 
characterized by: (1) thermal and structural finite elements 

formulated with common geometric discretization for full compatibility 
during the coupling of the analyses, (2) accurate nodal and element 
temperatures provided by improved thermal finite elements, and (3) 
accurate thermal loads for the structural finite element analysis to 
further improve accuracy of the structural response. 

Basic concepts and procedures of the integrated thermal-structural 
finite element analysis are described. New thermal finite elements 
for improved thermal analysis accuracy are developed. Thermal finite 
elements which yield exact nodal and element temperatures for one- 
dimensional linear steady-state heat transfer problems are presented. 
These thermal finite elements are formulated based upon using closed- 
form solutions of the governing differential equations. For general- 
heat transfer problems where closed-form solutions are not available. 
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improved thermal finite elements are developed by employing the 
nodeless variable formulation. The nodeless variable finite element 
uses extra unknowns as element variables to permit higher order 
element temperature interpolation functions. Detailed element 
temperature distributions are obtained without using additional 
element nodes while a common discretization with lower order congruent 
structural finite elements are maintained. 

Nodeless variable finite elements are formulated for the 
following heat transfer cases: (1) one-dimensional linear transient 

analysis, (2) one-dimensional nonlinear transient analysis with 
radiation, and (3) two-dimensional linear transient analysis. 

General formulations of the nodeless variable finite elements for 
each heat transfer case are described in detail. For comparison, 
conventional finite elements customarily used in typical finite 
element programs are also presented. Results of temperatures obtained 
from the thermal analysis are transferred directly to the structural 
analysis to compute displacements and stresses. 

To demonstrate the capabilities and efficiency of the integrated 
finite element approach, several examples in academic and more 
realistic problems are employed. The accuracy of the approach is 
evaluated by comparisons with analytical solutions and conventional 
thermal-structural analyses. Results indicate that the integrated 
finite element approach provides a significant improvement in the 
accuracy and efficiency of thermal-structural analysis and offers 
potential for applications to other coupled problems. 
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APPENDIX A 

EXACT FINITE ELEMENT INTERPOLATION FUNCTIONS 


Exact element interpolation functions in the form of equation 
(3.14) for the thermal finite element Cases 1-8 (Figure 4 and Table 1, 
pp. 39 and 4Q) are presented. Nodeless parameters are shown in 
Table 2 (p. 41)- The lower case letters in parentheses denote heat 
load cases defined in Table 1. General solutions to the differential 
equations for Cases 6 and 7 appear in reference [37]. 
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Hollow Cylinder (Case 3) 
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Conical Shell (Case 6) 
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Spherical Shell (Case 7) 
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ln l+sin(L/a) 
l-sin(L/a) 


(a,c,d) 


N q = In [cos (s/a)] - N 2 


ln [cos(L/a)] 
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Flow Passage (Case 8) 


N 


1 




(a,d) 



„ _ _ax sinh B(L-x) 

N 1 = 6 sinh 6L 


.. a(x-L) sinh 6x 

n s e 

w 2 sinh 0L 


(b) 


1 - - N 2 


where a = mc/2kA, 0 



+ m 


and m 


/hp/kA . 
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APPENDIX B 

FINITE ELEMENT MATRICES FOR ONE -DIMENSIONAL 
LINEAR STEADY-STATE PROBLEMS 

Exact finite element matrices for the thermal and structural 
finite elements described in Chapter 3 are presented. Thermal 
conductance matrices and heat load vectors are given in the form of 
equation (3.27) for Cases 1-6 and equation (3.17) for Case 8 
(Figure 4 and Table 1, pp. 39 and 40 ) . These finite element 
matrices are derived using the exact element interpolation functions 
shown in Appendix A. Similarly, structural stiffness matrices and 
equivalent nodal forces due to thermal loads are derived using the 
element displacement interpolation functions shown in Tables 3 and 
4 (pp. 61 and 67). The lower case letters in parentheses denote 
heat load cases defined in Table 1. 
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THERMAL FINITE ELEMENT MATRICES 

Rod (Case 1) 

Conductance Matrices: 



- K 22 - kA/L 

(a, c,d) 

*!2 

- _ kA/L 


*11 

” ^22 = cosh mL)/(m sinh mL) 

(b) 

*12 

= - hp/(m sinh mL) 


Eeat Load Vectors: 



Q x - Q 2 = QAL/2 

(c) 


Q x = Q 2 - qpL/2 

(d) 


where m = /hp/kA . 

Slab (Case 2) 

Conductance Matrices: 

Ku = K 22 = k/L (a,c) 

K 12 = - k/L 
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Heat Load Vector 

Q x = Q 2 = QL/2 


Hollow Cylinder (Case 3) 


Conductance Matrices 


Kii = K 22 = k/w 


^2 


= - k/w 


Heat Load Vector 



where w = ln(b/a) 


(- a 2 / 2 + (b 2 - a 2 ) / 4w) 
(b 2 /2 - (b 2 - a 2 ) / 4w ) 


Hollow Sphere (Case 4) 

Conductance Matrices 

Ku = K 22 = kab/(b-a) 

= - kab/ (b-a) 

Heat Load Vector 

Q l = Qa(2a 3 + b 3 - 3a 2 b) /(6(b-a)) 


(c) 


(a,c) 


(c) 


(a,c) 


(c) 
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Q 2 = Qb (a 3 + 2b 3 - 3ab 2 ) /(6(b-a)) 


Cylindrical Shell (Case 5) 


Conductance Matrices 


*11 = K 22 . k/L 


(a,c,d) 


K 12 = " k/L 


= K 22 = (h cosh mL) / (mt sinh mL) 


(b) 


K^2 - ~ h / (mt sinh mL) 


Heat Load Vectors 


Q x “ q 2 - QL/2 


Ql = q 2 = qL/2t 


= Q 2 = hT^ (cosh mL - 1) / (mt sinh mL) 


(c) 

(d) 

(b) 


where m - /h/kt . 


Conical Shell (Case 6) 


Conductance Matrices 


K n = K 22 - k/w 


(a,c,d) 


K^ 2 “ - k/w 
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Flow 


Heat Load Vectors 

= Q(-a 2 /2 + (b 2 -a 2 )/4w) (c) 

Q 2 = Q(b 2 /2 - (b 2 -a 2 )/4w) 

\ = q(-a 2 /2t + (b 2 -a 2 )/4wt) (d) 

Q 2 = a(b 2 /2t - (b 2 -a 2 )/4wt) 
where w = ln(b/a) 


Passage (Case 8) 


Conductance Matrices 


K c — K = kaa(e 2 “ L + l)/(e 2 “ L -1) (a,d) 

11 c 22 


K„ 


'12 


K„ 


•11 


K, 


'22 




'12 


= kA( (3H/2G) - (a/2) - (g 2 E(E+F) /2aG 2 ) ) 

= kA( (gH/2G) + (a/2) - (g 2 E(E-F) /2aG 2 ) ) 

= K c = kA(-(g 2 EH/2aG 2 ) - (3F/2G)) 

21 


=1^ = - mc/2 

11 12 


21 


K = mc/2 
V 22 


(b) 


(a,d) 
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K„ 


11 


12 


*h 


11 


Ku 


22 




12 


= - K = - mc/2 
22 


= - K__ = BE/ctG 

21 


= hp((BH/G) + (a) - (3 2 E(F+E)/aG 2 )) / 2 (g 2 -a 2 ) 


= hp((0H/G) - (a) - (B 2 E(F-E)/aG 2 )) / 2(g 2 -a 2 ) 


= K h = hp((B 2 EF/aG 2 ) - (BF/G)) / 2 (g 2 -a 2 ) 


(b) 


(b) 


Heat Load Vectors 

q i = qp (1 -e 2aL + 2aL) /2a(l -e 2aL ) (d) 

Q 2 = 3P(“1 +e -2aL e ) /2a(l -e ) 

Q 1 = hpT^ ( 3 (H+E-F ) -aG) / G(g 2 -a 2 ) (b) 

Q 0 = hpT (B(H-E-F) + aG) / G(B 2 -a 2 ) 

Z 00 

r 2 2 ~ 

where a = mc/2kA, 6 = + m , m = /hp/kA , E = sinh aL , 

F = cosh aL , G = sinh gL , H = cosh BL . 



175 


STRUCTURAL FINITE ELEMENT MATRICES 

T russ Element (Case 1) 

Stiffness Matrices 

K 11 = K 22 * (a,b,c,d) 

K 12 = K 22 ® 

Force Vectors 

F 1 = ~ F 2 = ~ a EA(T Q /6 + (a,c,d) 

Pi - - F 2 - - aEA(C 1 T 0 + C 2 (T l +T 2 )) (b) 

where C*j_ = 1 - (2 (cosh mL - 1) /mL sinh mL) 

= (cosh mL - l)/mL sinh mL 

Axisymmetric Element (Case 3, Plane Stress) 

Stiffness Matrices 

K 11 = E((b 2 +a 2 ) - v(b 2 - a 2 ))/(l -v 2 )(b 2 -a 2 ) (a,c) 

K 22 = E ^ 1>2 ~ a2 ^ + v (b 2 - a 2 )) /(I -v 2 )(b 2 - a!") 

K 12 = E(-2ab)/(l - v 2 )(b 2 -a 2 ) 
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Force 


where 


Vectors 

F 1 = - aKa P/2w(l - v) (b 2 - a' 
F 2 = bEaP/2w(l -V)(b 2 -a 2 ) 

P = ( - 2a 2 v + b 2 - a 2 ) (T x + a 2 w T / b 2 
+ (2b 2 w - b 2 + a 2 )(T 2 + w T q ) 

- (b 4 - a 4 ) w 2 T Q / b 2 - 2(b 2 - a 2 


) (a,c) 


) 


w T re f 
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APPENDIX C 

FINITE ELEMENT MATRICES FOR ONE-DIMENSIONAL 
LINEAR TRANSIENT PROBLEMS 


Finite element capacitance matrices for the thermal rod and 
axisymmetric elements described in Chapter 4 are presented. The 
conductance matrix coefficients Kqq , and heat load vector 
components are presented; the coefficients and Q^, 

i,j = 1,2 appear in Appendix B, Cases 1 and 3. The lower case 
letters in the parentheses denote heat load cases defined in Table 1. 
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Rod Element 

Capacitance Matrices 
C 00 = pcAL/30 

Cq^ = C ^2 = pcAL/12 

C 11 = C 22 = P cAL / 3 

C 12 = 

Cqq = PcA(((cosh mL - l)/sinh mL) (L/sinh mL - 3/m) + L) 

o 

Cqi = Cq 2 * pcA( (1 - cosh mL) (mL - sinh mL)/2m sinh mL) 

2 

C u = C 22 = pcA((sinh mL cosh mL - mL) /2m sinh mL) 

2 

C 12 = P 0 ^^ 11 ^ costl mL - sinh mL)/2m sinh mL) 

Conductance Matrices 

K 00 " kA/3L ( 

Kqq = (hp/m)(mL - 2 (cosh mL - l)/sinh mL) 

Heat Load Vectors 

Q q = QAL/6 


a,c,d) 


(b) 


a,c,d) 

(b) 

(c) 


Q 0 = qpL/6 


(d) 
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Qq - hpT^ (L - 2 (cosh mL - l)/m sinh raL) (b) 

where m = /hp/kA 

Axisymmetric Element 

Capacitance Matrices 

C Q0 = pc(b 2 -a 2 ) (4w 2 (a 4 + a 2 b 2 + b 4 ) + 9w(a 4 -b 4 ) (a,c) 

+ 6(a 2 -b 2 ) 2 )/24b 4 

C 01 = " P c < 4a4w2 + w(7a 2 +3b 2 ) (a 2 -b 2 ) + 4 (a 2 -b 2 ) 2 ) /16wb 2 
C Q2 = pc (4b 4 w 2 + w(7b 2 +3a 2 ) (a 2 -b 2 ) + 4(a 2 -b 2 ) 2 ) /16wb 2 

C n = pc (b 2 -a 2 (l+2w+2w 2 ) ) /4w 2 

2 2 2 2 2 
C 12 = pc(a -b + w(a -b ))/4w 

C 22 = pc(b 2 (l-2w+2w 2 ) - a 2 )/4w 2 

Conductance Matrices 

K Q0 = kw<w(l-(a/b) 4 ) - (l-(a/b) 2 ) 2 ) (a,c) 

Heat Load Vector 

Q q = (Qb 2 /4)(w(l-a 4 /b 4 ) - (l-a 2 /b 2 ) 2 ) 


(c) 
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APPENDIX D 

FINITE ELEMENT MATRICES FOR ONE -DIMENSIONAL NONLINEAR 
TRANSIENT ANALYSIS WITH RADIATION HEAT TRANSFER 

Jacobian matrices and residual heat load vectors for the 
conventional finite element and the nodeless variable finite 
element described in Chapter 5 are presented. These element 
matrices which appear in Eq. (5.34) are given in the form of computer 
subroutines. The subroutines are written in FORTRAN IV where the 
definitions of variables used are provided. 



r>r>or><->or>or>oor>or>r>or>r>or>oor>r>or>or>r>r>oooooo 


**************************** 

VARIABLES USED IN THE FOLLOWING SUBROUTINES ARE 
AS FOLLOUSi 


* * * * 
DEFINED 


TK 

MATERIAL THERMAL CONDUCTIVITY 




RHO 

DENSITY 





CP 

SPECIFIC HEAT 





AREA 

ROD CROSS-SECTIONAL AREA 





XL 

ELEMENT LENGTH 





PS 

CROSS-SECTIONAL PERIMETER 

FOR 

EMITTED 

ENERGY 

PQ 

CROSS-SECTIONAL PERIMETER 

FOR 

INCIDENT 

ENERGY 

EMIS 

SURFACE EMISSIVITY 





ABSORP 

SURFACE ABSORPTIVITY 





STEFAN 

STEFAN-BOLTZMANN CONSTANT 





QDOT 

INCIDENT HEATING RATE PER 

UNIT 

AREA 



TM ( ) 

ELEMENT NODAL TEMPERATURES 

AT 

THE M 

ITERATION 

TN ( ) 

ELEMENT NODAL TEMPERATURES 

AT 

THE N 

TIME STEP 

DELTA 

TIME INCREMENT USED IN CRANK-NICOLSON 

ALGORITHM 

ICONS 

• EQ • 1 CONSISTENT FORMULATION 





• NE • 1 LUMPED FORMULATION 


ALL JACOBIAN MATRICES ARE REPRESENTED BY VARIABLES BEGIN WITH 

AJ (_#_) E.G. A JR AD REPRESENTS JACOBIAN MATRIX CONTRIBUTED 

FROM CONDUCTANCE RADIATION MATRIX. 

ALL RESIDUAL LOAD VECTORS ARE REPRESENTED BY VARIABLES BEGIN 

WITH R <_) E.G. RR AD REPRESENTS RESIDUAL HEAT LOAD VECTOR 

FROM CONDUCTANCE RADIATION MATRIX. 
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c ******************************** 

c 

c 

SUBROUTINE JRCOND < TK# ARE A# XL , TM, AJCOND, RCONO ) 

C 

C SET UP JACOBIAN CONDUCTION MATRIX AND CONDUCTION 

C RESIDUAL LOAD VECTOR FOR CONVENTIONAL ROD ELEMENT 

C 

DIMENSION AJCOND(2»2),RCOND(2),TM(2) 

C 

XX ■ TK+APEAf XL 
AJCOND (1,1) • XX 
AJCOND (1*2) ■ -XX 
AJCOND ( 2 f 1 ) • -XX 
A JCONDt 2# 2) - XX 
C 

RCOND(l) • -XXMTM(1)-TM<2)) 

RCOND ( 2 ) ■ - RCOND(l) 

C 

RETURN 

END 

C 

C 

C ******************************** 

C 

C 

SUBROUTINE JRR AO ( EMI S , ST EF AN, PS , XL > TM, ICONS t A J R AD» RR AD ) 

C 

C SET UP RADIATION JACOBIAN MATRIX AND RADIATION RESIDUAL 

C LOAD VECTOR FOR CONVENTIONAL ROD ELEMENT 

C 

DIMENSION AJRAD(2»2)»RRAD(2)#TM(2 ) 

C 

IF ( ICONS ♦ NE • 1 ) 60 TO 10 

XX • EMIS*STEFAN*PS*XL/15. 

T1 • TM ( 1 ) 
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T2 • TM ( 2 ) 

T1S » T1*T1 
T2S • T2*T2 
TIC • T1S+T1 
T2C ■ T2S+T2 
C 

T1ST2 • T1S+T2 
T1T2S ■ T1*T2S 

AJRAO(lfl) ■ XX* ( 10»*T1C t 6.*T1ST2 ♦ 3.*T1T2S t T2C ) 

A J R AD ( 1 * 2 ) ■ XX* f 2«*T1C * 3.*T1ST2 ♦ 3.*T1T2S ♦ 2.*T2C) 

A JR AD ( 2» 1 ) - A JR AD ( 1 » 2 ) 

A JRAD(2»2 ) ■ X X* ( TIC * 3.*T1ST2 ♦ 6.+T1T2S ♦ 10.*T2C> 

C 

T1F • T1C*T1 
T2F • T2C *T2 

C 

T1CT2 % T1C*T2 
T1ST2S- T1S+T2S 
T1T2C % T1*T2C 
YY % —XX/ 2 • 

RR AD ( 1 ) % YY*(5.*T1F ♦ <t.*TlCT2 ♦ 3.*T1$T2S + 2. + T1T2C ♦ 
RRADC2) ■ YY* ( T1F ♦ 2.*T1CT2 ♦ 3.+T1ST2S ♦ A.+T1T2C ♦ 5 

C 

RETURN 

C 

10 CONTINUF 

XX ■ 2.*FMIS*STEFAN*PS*Xl 
A JR AD ( If 1 ) ■ XX* ( TM ( 1 ) ♦* 3 • ) 

AJR AOI 1# 2 l - 0. 

A JR AD( 2* 1 ) % 0. 

A JR AD( 2»2 ) ■ XX*(TM<2>**3.) 

C 

YY • -XX/A. 

RRAD(l) ■ YY*(TM(l)* + <».) 

RR AD ( 2 ) • YY*(TM(2)**4.) 

C 


T2F ) 
• *T2F ) 
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r> o r> r> r> r> o o o <r> o r> r> o r> o r> r> o o 


RETURN 

END 

4 + 4 + + + + + + + + *** + + + ** + + + *** + + + + ** + 


SUBROUTINE RQDOT ( QDOT# AB SOPP# PQ» XL # RQ) 

SET UP INCIDENT RESIDUAL LOAD VECTOR FOR 
CONVENTIONAL ROD ELEMENT 

DIMENSION RQ I 2 ) 

XX • QDOT+ABSOPPARQ+XL/2. 

RQ(1) « XX 
RQ { 2 ) - XX 

RETURN 

END 


+ + + + + + + + + + + + + + + + + + + * + + + + + + + ++ + + + 


SUBROUTINE JRC AP < RHO# C P# AR E A# XL # DE LT A# TM# TN# ICONS# A JCAP#RCAP ) 

SET UP CAPACITANCE JACOBIAN MATRIX AND CAPACITANCE RESIDUAL 
LOAD VECTOR FOR CONVENTIONAL ROD ELEMENT 

DIMENSION AJCAP(2#2)#RCAP(2)#TM(2) #TN<2) 

IF(ICONS.NE.l) 60 TO 10 

XX - RH0*CP+AREA*XL/(3«*DELTA) 

AJCAP<1»1) ■ 2 • + XX 
A JC AP ( 1# 2 ) - XX 


184 


r> o n o r> r> o r> 


AJCAP(2#1) ■ A JC AP ( 1# 2 ) 

AJCAP(2,2) - AJCAP(1>1) 

C 

RCAP(l) - XX*(2.+TN<1) - 2 . *TH « 1 > ♦ TN(2) - TM(2)» 

RCAP (2 ) ■ XX*( TN ( 1 ) - TH < 1 1 + 2.*TN(2) - 2.*TH(2)> 

C 

RETURN 

C 

10 CONTINUE 

YY - RHO*CP*AREA* XL /DELTA 
AJCAPI1,1) ■ YY 
A JC AP ( 1» 2 ) • 0. 

A JC AP ( 2> 1 ) • 0. 

A JC AP ( 2# 2 ) ■ YY 
C 

RCAP(l) ■ YY+( TN ( 1 )-TM(l 1 1 
RC AP (2 ) • YY* ( TN < 2 l-TH (2 ) I 
C 

PETURN 

END 




SUBROUTINE I JRCOND < TK , AR E A# XL# TM# A JC OND* RC OND ) 

SET UP JACOBIAN CONDUCTION MATRIX AND CONDUCTION 
RESIDUAL LOAD VECTOR FOR NODELESS VARIABLE ROD ELEMENT 

DIMENSION AJC0ND(3«3)»RC0ND(3)»TM(3) 

C 

XX ■ TK+AREA/XL 
AJCOND ( 1# 1 ) • XX 
AJC0ND(1,2) • -XX 
AJCOND ( 1# 3 ) • 0. 
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o r> r> 


AJC0ND<2,1) • -XX 
A JCOND (2,2) - XX 
A JCOND( 2,3 ) • 0. 

A JCOND ( 3» 1 ) ■ 0. 

A JCOND (3,2) - 0. 

A JCOND ( 3# 3 ) - X X / 3 • 

C 

PCONDU) • -XX* ( TM ( 1 )-TM ( 2 ) ) 

RCOND( 2 ) • -RCONO(l) 

RCQNO( 3 ) • -XX*TH(3)/3. 

C 

RETURN 

END. 

C 

C 

C ******************************** 

C 

C 

SUBROUTINE I JRR AD ( EM I S ,S TEF AN, PS , X L, TM, A JR AD, RP AD ) 

C 

SET UP RADIATION JACOBIAN MATRIX AND RADIATION RESIDUAL 
LOAD VECTOR FOR NODELESS VARIABLE ROD ELEMENT 

DIMENSION AJRAD(3,3)»RRAD(3),TM(3) 

C 

XX ■ EMIS*STEFAN»PS*XL/630. 

T 1 ■ TM ( 1 ) 

T2 ■ TM ( 2 ) 

TO - TM ( 3 ) 

TIS ■ T1*T1 
T2S • T2*T2 
TOS ■ TO+TO 
TIC • T1S+T1 
T2C - T2S+T2 
TOC - TOS ♦TO 
C 
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c 


1 

1 

1 

1 

1 

1 


T0T1 - T0*T1 
TOT1S - T0*T1S 
TOST1 ■ T0S*T1 

A JR AO ( If 1 ) • XX*(42.*T2C + (126. *T1 ♦ 54.*TO)*T2S 

+ (252.*T1S + 144 .*TOTl + 27.*TOS)*T2 

♦ 420.*T1C ♦ 180.+T0T1S + 45.*T0ST1 
+ 5,*T0C 

A JR AD ( 1 f 2 ) ■ XX* ( 84.*T2C ♦ (126. *T1 ♦ 72.*T0)*T2S 

♦ (126. + T1S t 108 ,*T0T1 ♦ 27.*T0S)*T2 

♦ 84 • *T 1C ♦ 72 »*T0T1S ♦ 27.*T0ST1 

♦ 4 , *TOC 

A JR AD ( If 3 ) - XX*(24.*T2C ♦ (54.*T1 + 27.*T0)*T2S 


1 

1 

1 

AJRAD(2»1) 

AJRAD(2#2) 

1 

1 

1 


♦ (72.+T1S ♦ 54.*T0T1 ♦ 12.*T0S)*T2 

♦ 60 • *T1C ♦ 45.*T0T1S ♦ 15.+T0ST1 

♦ 2. *TOC 
A JR AD ( If 2 ) 

XX* ( 420,*T2C * (252. *T1 ♦ 180.*T0)*T2S 

♦ ( 126. *T1S ♦ 144 **T0T1 ♦ 45.*TOS)*T2 

♦ 42.*T1C ♦ 54 ,*T0T1S + 27.*T0ST1 

♦ 5.*T0C 


AJR AD( 2# 3 ) 

1 

1 

1 

AJRAD(Bfl) 
AJRADOf 2) 
A JR AD ( 3# 3 ) 


1 

1 

1 


XX* ( 60. + T2C ♦ ( 7 2 • * T 1 ♦ 45.*T0)*T2S 

+ ( 54 .*T1S ♦ 54.*T0T1 ♦ 15.*T0S)*T2 
+ 24.+T1C + 27»*T0T1S + 12.*T0ST1 

♦ 2 , *TOC 
A JR AD( If 3 ) 

A JR AD ( 2f 3 ) 

XX* ( 15.*T2C ♦ ( 27 . *T1 ♦ 15.*T0)*T2S 

♦ ( 27 .*T1S ♦ 24 .* T0T1 ♦ 6.*T0S)*T2 

♦ 15 ,*T1C ♦ 15.*T0T1S ♦ 6 . *T0ST1 
+ 10.*T0C/11. 


) 


) 


) 


) 


) 


) 


T1F - TIC *T1 
T2F - T2C *T2 
TOF - T OC * TO 
T0T1C - T0*T1C 
T0ST1S ■ T0S+T1S 
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r> o r> r> o 


TOCTl • TOC*Tl 

C 

YY ■ -XX/2. 

RR AD ( 1 ) • YY*(42.*T2F ♦ 184. *T1 ♦ 48.*T0)*T2C 
1 ♦ ( 1 26. *T1S + 108.+T0T1 ♦ 27.*TOS)*T2S 

1 ♦ ( 168 • +T1C + 144.+T0T1S + 54.*TOSTl ♦ 8.*T0C)*T2 

1 + 210.*T1F + 120.+T0T1C ♦ 45.*T0ST1S + 10.+T0CT1 ♦ TOF ) 

RR AD ( 2 ) ■ YY* (210. *T2F ♦ (168. +T1 ♦ 120,*T0)+T2C 
1 ♦ 1126. *T1S + 144.+T0T1 ♦ 45.*T0$I*T2S 

1 ♦ (B4.+T1C ♦ 108 . *T0T1 S ♦ 54.+T0ST1 ♦ 10.*T0C)*T2 

1 ♦ 42.*T1F ♦ 4B.*T0T1C ♦ 27.*T0ST1S ♦ 8.+T0CT1 ♦ TOF ) 

RR AD (3 1 • YY*(30.*T2F ♦ (40.+T1 ♦ 30.*T0)*T2C 
1 ♦ ( 5 4. * T1 S ♦ 54. + T0T1 ♦ 15.*T0S)*T2S 
1 + I48.+T1C ♦ 54.+T0T1S ♦ 24.*T0ST1 ♦ 4.*T0C)*T2 

1 ♦ 30 • *T IF ♦ 30. *T0T1C ♦ 15.*T0ST1S ♦ 4.+T0CT1 ♦ 5.*T0F/11.) 

C 

RETURN 

END 


******************************** 


SUBROUTINE IRQDOT ( ODOT #ABSORP»PQ»XL#PQI 
C 

C SET UP INCIDENT RESIDUAL LOAD VECTOR FOR 

C NODELESS VARIABLE ROD ELEMENT 

C 

DIMENSION RQ ( 3) 

C 

XX ■ QDOT*ABSORP*PO*XL 
R Q ( 1 ) - XX/2. 

PQ ( 2 ) - R 0 < 1 ) 

RO ( 3 ) • X X / 6 • 

C 

RETURN 
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r> o o r> o o r> o o o r> o r» o r> 


END 


****************************4*** 


SUBROUTINE I JRC AP ( RHOfCP f AREA# XL fDELTAf TH, TNf A JC APf PC AP ) 

SET UP CAPACITANCE JACOBIAN MATRIX AND CAPACITANCE RESIDUAL 
LQAO VECTOR FOR NODELESS VARIABLE ROD ELEMENT 


DIMENSION AJCAPI3f3)fRCAPI3)fTMI3»f TNI3I 


XX • 2**RHO*CP+AREA*XL/DELTA 


AJC AP ( 1 * 1 ) 
A JCAPC 1# 2 ) 
A JC AP ( I# 3 ) 
AJCAP|2#1) 
A JC AP ( 2* 2 ) 
AJC API 2#3 ) 
A JC AP 1 3f 1 ) 
A JC AP I 3# 2 ) 
AJC AP ( 3# 3 ) 


X X / 3 • 

XX/6. 

XX/ 12. 
AJCAPIlf 2) 
A JC AP « 1# 1 > 
AJCAPIl# 3) 
AJCAPU, 3) 
AJCAP(2f 3) 
XX/30. 


DTO - TNI 3 »-TM 13) 

DTI • TNI 1 l-TMIl ) 

DT2 - TN I 2 » -TH I 2 1 

RCAP(l) - XX+IDTO/12. ♦ DTI/3. + 0T2/6. ) 

RCAPI2) - XX+IDTO/12. ♦ DTI/6. ♦ DT2/3. ) 

RCAPI3) ■ XX* I DT0/30* ♦ DTI/12. ♦ DT2/12.I 


RETURN 

END 


*********** 


********************* 
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